Preliminary NPB Model for Birth Weight
Sheena Martenies
06/26/2020
1 Exploring the data set
The HS data set was previously used in the CEI paper (Martenies et al., 2019). In the original analysis, we used an exposure index based on the CalEnvironScreen tool. We observed lower birth weights and lower adiposity associated with higher index scores, driven largely by exposures to social indicators of health at the neighborhood level. Now, we are aiming to use methods for mixtures to try to identify which exposures are driving these association.
The complete data set for the birth weight outcome consists of n = 897 participants. This represents 63.62% of the original Healthy Start 1 cohort.
Of the 897 participants, 0.27% identify as Latina, 0.17% identify as Black, and 0.27% identify as another non-NHW race or ethnicity. The median age of mothers in this dataset is 28 years. 0.51% of babies born were male.
1.1 Exposure data
We have included 19 exposures in our analysis.
These exposures are based on the census tract where each mother lived at the time of enrollment into Healthy Start. With the exception of air pollution (mean_pm and mean_o3), these are based on long-term averages at for each census tract. For mean_pm and mean_o3 are based on the average pollution levels across each pregnancy (est. conception date to delivery date) estimated using ordinary kriging and monitoring data.
#' Exposure data
X <- select(hs_data2, mean_pm, mean_o3, pct_tree_cover, pct_impervious,
mean_aadt_intensity, dist_m_tri:dist_m_mine_well,
cvd_rate_adj, res_rate_adj, violent_crime_rate, property_crime_rate,
pct_less_hs, pct_unemp, pct_limited_eng, pct_hh_pov) %>%
as.matrix()
head(X)## mean_pm mean_o3 pct_tree_cover pct_impervious mean_aadt_intensity
## [1,] 8.483046 47.19072 6.006276 43.30893 10128.4962
## [2,] 6.598608 50.05090 7.281109 48.36432 10749.0359
## [3,] 7.454146 48.57052 17.205991 31.67281 9048.6468
## [4,] 6.671239 50.06429 6.842898 45.00359 4223.3434
## [5,] 7.122537 50.14275 3.357792 28.16745 858.7283
## [6,] 7.637453 47.03125 10.743612 45.87564 15603.9800
## dist_m_tri dist_m_npl dist_m_waste_site dist_m_major_emit dist_m_cafo
## [1,] 2827.538 729.2371 4829.780 7968.654 29116.58
## [2,] 1576.420 5239.2211 4417.792 3780.951 51044.30
## [3,] 3350.303 2992.2968 5211.871 7423.232 36079.21
## [4,] 3364.954 6998.1286 8921.318 9636.816 42235.78
## [5,] 2923.811 3427.2247 7006.042 6806.912 29145.98
## [6,] 3364.200 3166.5395 4484.960 5265.285 43921.85
## dist_m_mine_well cvd_rate_adj res_rate_adj violent_crime_rate
## [1,] 1749.1256 275.2480 155.7767 14.377133
## [2,] 7354.5310 279.6435 226.8038 8.905404
## [3,] 4887.2996 221.0414 157.6974 7.636888
## [4,] 3752.6399 203.8812 142.5368 2.850212
## [5,] 729.7784 194.1983 101.0046 5.435988
## [6,] 5870.6867 174.3361 120.3281 5.035971
## property_crime_rate pct_less_hs pct_unemp pct_limited_eng pct_hh_pov
## [1,] 37.32935 31.784946 11.529628 26.114650 12.010919
## [2,] 67.03932 15.290231 4.908306 8.500401 18.123496
## [3,] 46.78194 6.891702 4.564963 0.000000 6.307978
## [4,] 21.95270 2.725915 5.623583 1.350621 9.292274
## [5,] 22.49834 12.919186 5.234103 6.307385 2.115768
## [6,] 47.15500 3.842365 10.000000 5.121799 25.171768Variance and histograms of the exposure variables (in their original units):
var(X)## mean_pm mean_o3 pct_tree_cover
## mean_pm 0.391784015 0.006083605 -0.2054297
## mean_o3 0.006083605 9.383489039 -0.4158151
## pct_tree_cover -0.205429726 -0.415815089 9.7193077
## pct_impervious 0.508898445 -1.674151031 5.8719893
## mean_aadt_intensity -182.234953786 474.627052967 8431.6446632
## dist_m_tri -255.176839682 444.286548683 -73.1423054
## dist_m_npl -289.002141382 539.849185829 434.4654007
## dist_m_waste_site -275.262105884 261.902915064 1933.8647304
## dist_m_major_emit 71.096593638 577.257325397 265.4284518
## dist_m_cafo -1291.237441927 -35.275020052 10170.6234275
## dist_m_mine_well -339.250592215 -375.434990683 3136.3680766
## cvd_rate_adj 3.871688575 0.939328342 -24.8232924
## res_rate_adj 2.356328835 -0.181515705 -3.5331376
## violent_crime_rate 0.232839920 0.577648302 -4.0583754
## property_crime_rate 2.001989749 -2.773092354 -22.6429724
## pct_less_hs 1.132232861 0.637361326 -7.5753471
## pct_unemp 0.100439902 0.288530482 -0.3330523
## pct_limited_eng 0.432516169 0.295617023 -2.8349116
## pct_hh_pov 0.731824476 -0.606648513 0.3805472
## pct_impervious mean_aadt_intensity dist_m_tri
## mean_pm 0.5088984 -182.2350 -255.17684
## mean_o3 -1.6741510 474.6271 444.28655
## pct_tree_cover 5.8719893 8431.6447 -73.14231
## pct_impervious 176.8316214 55459.6063 -15279.44024
## mean_aadt_intensity 55459.6063235 67283287.0201 -1315386.69307
## dist_m_tri -15279.4402428 -1315386.6931 6558190.20296
## dist_m_npl -7729.3843793 1683196.0799 4282727.94125
## dist_m_waste_site -4662.9983638 2039577.9230 2441267.84540
## dist_m_major_emit 2627.0270993 2477155.3406 1433153.16531
## dist_m_cafo 16586.9964129 15462371.9832 3431065.70215
## dist_m_mine_well 706.6674650 2073244.5987 995872.11873
## cvd_rate_adj 230.4542985 20477.4374 -49347.60273
## res_rate_adj 176.8108084 33055.3733 -31870.98664
## violent_crime_rate 26.6945028 5736.5627 -1014.08753
## property_crime_rate 118.0737725 22077.3894 -5365.69285
## pct_less_hs 56.8383947 -4056.6889 -12372.14262
## pct_unemp 25.9434246 6003.3343 -2527.22451
## pct_limited_eng 41.9919053 2620.6198 -5408.86434
## pct_hh_pov 82.2198624 17850.1649 -8842.76408
## dist_m_npl dist_m_waste_site dist_m_major_emit
## mean_pm -289.0021 -275.2621 71.09659
## mean_o3 539.8492 261.9029 577.25733
## pct_tree_cover 434.4654 1933.8647 265.42845
## pct_impervious -7729.3844 -4662.9984 2627.02710
## mean_aadt_intensity 1683196.0799 2039577.9230 2477155.34057
## dist_m_tri 4282727.9413 2441267.8454 1433153.16531
## dist_m_npl 11125411.7474 4193498.0586 6948817.25739
## dist_m_waste_site 4193498.0586 5344101.7540 1395277.06805
## dist_m_major_emit 6948817.2574 1395277.0681 10114549.72263
## dist_m_cafo 5416531.1320 5586018.8251 -2993791.05377
## dist_m_mine_well 256924.3029 1375784.7856 -1810174.74785
## cvd_rate_adj -30921.0390 -43119.5785 16272.40152
## res_rate_adj -19393.1304 -32402.8440 -1320.21297
## violent_crime_rate -672.9264 -3702.6112 -360.49700
## property_crime_rate -18283.4264 -22350.3006 -24007.42305
## pct_less_hs -6760.5337 -11422.4985 8866.74917
## pct_unemp 2195.0515 -1476.4094 5212.74830
## pct_limited_eng 498.0033 -4277.8134 9367.28435
## pct_hh_pov -1135.3843 -7599.7432 8682.26135
## dist_m_cafo dist_m_mine_well cvd_rate_adj
## mean_pm -1291.23744 -339.2506 3.8716886
## mean_o3 -35.27502 -375.4350 0.9393283
## pct_tree_cover 10170.62343 3136.3681 -24.8232924
## pct_impervious 16586.99641 706.6675 230.4542985
## mean_aadt_intensity 15462371.98316 2073244.5987 20477.4373759
## dist_m_tri 3431065.70215 995872.1187 -49347.6027339
## dist_m_npl 5416531.13199 256924.3029 -30921.0389720
## dist_m_waste_site 5586018.82514 1375784.7856 -43119.5785165
## dist_m_major_emit -2993791.05377 -1810174.7478 16272.4015197
## dist_m_cafo 46324000.89481 9345575.3772 -46645.9665229
## dist_m_mine_well 9345575.37722 4430024.9964 -39046.5984701
## cvd_rate_adj -46645.96652 -39046.5985 2039.8569530
## res_rate_adj -13772.40263 -16322.5110 1289.5661935
## violent_crime_rate 722.31907 -2032.3464 135.9487143
## property_crime_rate -15833.92381 -4272.3829 343.9364726
## pct_less_hs -26060.83378 -10037.6577 328.3044447
## pct_unemp -1030.96916 -2827.2369 105.0153846
## pct_limited_eng -7089.15821 -4814.6687 183.5853966
## pct_hh_pov -855.38016 -5030.4055 266.1004715
## res_rate_adj violent_crime_rate property_crime_rate
## mean_pm 2.3563288 0.2328399 2.001990
## mean_o3 -0.1815157 0.5776483 -2.773092
## pct_tree_cover -3.5331376 -4.0583754 -22.642972
## pct_impervious 176.8108084 26.6945028 118.073773
## mean_aadt_intensity 33055.3733277 5736.5627383 22077.389365
## dist_m_tri -31870.9866403 -1014.0875345 -5365.692846
## dist_m_npl -19393.1304345 -672.9263612 -18283.426420
## dist_m_waste_site -32402.8439544 -3702.6111771 -22350.300554
## dist_m_major_emit -1320.2129699 -360.4970006 -24007.423046
## dist_m_cafo -13772.4026269 722.3190727 -15833.923813
## dist_m_mine_well -16322.5110008 -2032.3464340 -4272.382880
## cvd_rate_adj 1289.5661935 135.9487143 343.936473
## res_rate_adj 1091.1856742 104.4979610 333.780710
## violent_crime_rate 104.4979610 40.1175363 160.725724
## property_crime_rate 333.7807097 160.7257236 1295.004010
## pct_less_hs 197.8827546 22.5579950 -3.138375
## pct_unemp 72.3576933 11.3130282 1.362247
## pct_limited_eng 104.0524036 12.7978322 -14.963510
## pct_hh_pov 201.6582659 29.1947400 64.236239
## pct_less_hs pct_unemp pct_limited_eng pct_hh_pov
## mean_pm 1.1322329 0.1004399 0.4325162 0.7318245
## mean_o3 0.6373613 0.2885305 0.2956170 -0.6066485
## pct_tree_cover -7.5753471 -0.3330523 -2.8349116 0.3805472
## pct_impervious 56.8383947 25.9434246 41.9919053 82.2198624
## mean_aadt_intensity -4056.6889048 6003.3343312 2620.6197529 17850.1649192
## dist_m_tri -12372.1426191 -2527.2245090 -5408.8643368 -8842.7640785
## dist_m_npl -6760.5337115 2195.0514738 498.0033420 -1135.3843390
## dist_m_waste_site -11422.4985495 -1476.4094188 -4277.8133935 -7599.7432386
## dist_m_major_emit 8866.7491706 5212.7483023 9367.2843472 8682.2613524
## dist_m_cafo -26060.8337755 -1030.9691591 -7089.1582114 -855.3801591
## dist_m_mine_well -10037.6576614 -2827.2368665 -4814.6687400 -5030.4055237
## cvd_rate_adj 328.3044447 105.0153846 183.5853966 266.1004715
## res_rate_adj 197.8827546 72.3576933 104.0524036 201.6582659
## violent_crime_rate 22.5579950 11.3130282 12.7978322 29.1947400
## property_crime_rate -3.1383751 1.3622468 -14.9635105 64.2362387
## pct_less_hs 162.1681017 39.4206217 85.1910014 100.9072175
## pct_unemp 39.4206217 24.6546969 25.2172769 36.9693212
## pct_limited_eng 85.1910014 25.2172769 68.6532943 67.2758215
## pct_hh_pov 100.9072175 36.9693212 67.2758215 119.7157808ggplot(pivot_longer(as.data.frame(X), mean_pm:pct_hh_pov, names_to = "exp", values_to = "value")) +
geom_histogram(aes(x = value)) +
facet_wrap(~ exp, scales = "free")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Scaling the exposure variables
X.scaled <- apply(X, 2, scale)
head(X.scaled)## mean_pm mean_o3 pct_tree_cover pct_impervious mean_aadt_intensity
## [1,] 1.60876944 -0.2502907 -0.08261827 0.2084897 -0.02626143
## [2,] -1.40186806 0.6834152 0.32629926 0.5886571 0.04938980
## [3,] -0.03503482 0.2001460 3.50981981 -0.6665506 -0.15790801
## [4,] -1.28583023 0.6877893 0.18573791 0.3359288 -0.74617032
## [5,] -0.56482237 0.7133998 -0.93215008 -0.9301550 -1.15635722
## [6,] 0.25782234 -0.3023500 1.43693690 0.4015075 0.64126566
## dist_m_tri dist_m_npl dist_m_waste_site dist_m_major_emit dist_m_cafo
## [1,] -0.4008664 -1.44212141 -0.172775897 -0.1090638 -1.1354079
## [2,] -0.8894134 -0.08999604 -0.350992130 -1.4258114 2.0863310
## [3,] -0.1967327 -0.76363997 -0.007492447 -0.2805619 -0.1124210
## [4,] -0.1910117 0.43733697 1.597126043 0.4154599 0.7921355
## [5,] -0.3632728 -0.63324550 0.768623456 -0.4743524 -1.1310891
## [6,] -0.1913062 -0.71140077 -0.321936759 -0.9590895 1.0398621
## dist_m_mine_well cvd_rate_adj res_rate_adj violent_crime_rate
## [1,] -0.7748959 0.6917151 -0.2847192 0.2444833
## [2,] 1.8883051 0.7890362 1.8654606 -0.6194047
## [3,] 0.7160914 -0.5084805 -0.2265752 -0.8196807
## [4,] 0.1769998 -0.8884275 -0.6855261 -1.5754111
## [5,] -1.2592010 -1.1028189 -1.9428184 -1.1671633
## [6,] 1.1833114 -1.5425892 -1.3578439 -1.2303189
## property_crime_rate pct_less_hs pct_unemp pct_limited_eng pct_hh_pov
## [1,] -0.5008789 1.20091488 0.36793456 2.15383825 -0.3020187
## [2,] 0.3247153 -0.09436047 -0.96557105 0.02798425 0.2566428
## [3,] -0.2382061 -0.75386917 -1.03471894 -0.99792447 -0.8232412
## [4,] -0.9281724 -1.08099463 -0.82151747 -0.83491873 -0.5504902
## [5,] -0.9130100 -0.28055077 -0.89995692 -0.23668961 -1.2063898
## [6,] -0.2278392 -0.99332353 0.05987409 -0.37977738 0.9008223Variance and histograms of the exposure variables (scaled):
var(X.scaled)## mean_pm mean_o3 pct_tree_cover pct_impervious
## mean_pm 1.000000000 0.003172893 -0.105274276 0.06114033
## mean_o3 0.003172893 1.000000000 -0.043541201 -0.04109912
## pct_tree_cover -0.105274276 -0.043541201 1.000000000 0.14164056
## pct_impervious 0.061140333 -0.041099117 0.141640558 1.00000000
## mean_aadt_intensity -0.035493982 0.018889337 0.329716765 0.50844411
## dist_m_tri -0.159193706 0.056635532 -0.009161339 -0.44867872
## dist_m_npl -0.138426634 0.052836279 0.041781063 -0.17426369
## dist_m_waste_site -0.190232937 0.036984589 0.268331135 -0.15168685
## dist_m_major_emit 0.035715124 0.059253499 0.026770503 0.06211712
## dist_m_cafo -0.303095664 -0.001691929 0.479321466 0.18326720
## dist_m_mine_well -0.257510042 -0.058230361 0.477976199 0.02524829
## cvd_rate_adj 0.136954783 0.006789463 -0.176295758 0.38371166
## res_rate_adj 0.113962827 -0.001793836 -0.034307854 0.40251263
## violent_crime_rate 0.058730941 0.029772424 -0.205526308 0.31693831
## property_crime_rate 0.088879773 -0.025156291 -0.201827442 0.24673907
## pct_less_hs 0.142046220 0.016338825 -0.190810449 0.33564418
## pct_unemp 0.032317153 0.018969666 -0.021515177 0.39291400
## pct_limited_eng 0.083396591 0.011647067 -0.109746623 0.38111402
## pct_hh_pov 0.106858205 -0.018100029 0.011156172 0.56509450
## mean_aadt_intensity dist_m_tri dist_m_npl
## mean_pm -0.03549398 -0.159193706 -0.13842663
## mean_o3 0.01888934 0.056635532 0.05283628
## pct_tree_cover 0.32971677 -0.009161339 0.04178106
## pct_impervious 0.50844411 -0.448678724 -0.17426369
## mean_aadt_intensity 1.00000000 -0.062619247 0.06152095
## dist_m_tri -0.06261925 1.000000000 0.50138396
## dist_m_npl 0.06152095 0.501383960 1.00000000
## dist_m_waste_site 0.10755964 0.412369055 0.54385239
## dist_m_major_emit 0.09495686 0.175965418 0.65505772
## dist_m_cafo 0.27696155 0.196849357 0.23859436
## dist_m_mine_well 0.12008641 0.184760224 0.03659688
## cvd_rate_adj 0.05527416 -0.426652406 -0.20525615
## res_rate_adj 0.12199419 -0.376750794 -0.17601130
## violent_crime_rate 0.11041575 -0.062519618 -0.03185242
## property_crime_rate 0.07479259 -0.058223492 -0.15232247
## pct_less_hs -0.03883608 -0.379376300 -0.15916230
## pct_unemp 0.14739715 -0.198747643 0.13253691
## pct_limited_eng 0.03855846 -0.254907958 0.01801953
## pct_hh_pov 0.19888999 -0.315587892 -0.03111065
## dist_m_waste_site dist_m_major_emit dist_m_cafo
## mean_pm -0.19023294 0.03571512 -0.303095664
## mean_o3 0.03698459 0.05925350 -0.001691929
## pct_tree_cover 0.26833114 0.02677050 0.479321466
## pct_impervious -0.15168685 0.06211712 0.183267205
## mean_aadt_intensity 0.10755964 0.09495686 0.276961553
## dist_m_tri 0.41236906 0.17596542 0.196849357
## dist_m_npl 0.54385239 0.65505772 0.238594356
## dist_m_waste_site 1.00000000 0.18977973 0.355027509
## dist_m_major_emit 0.18977973 1.00000000 -0.138307324
## dist_m_cafo 0.35502751 -0.13830732 1.000000000
## dist_m_mine_well 0.28275484 -0.27042332 0.652378929
## cvd_rate_adj -0.41298787 0.11328659 -0.151743889
## res_rate_adj -0.42432287 -0.01256670 -0.061257230
## violent_crime_rate -0.25287368 -0.01789622 0.016755559
## property_crime_rate -0.26866460 -0.20976678 -0.064647236
## pct_less_hs -0.38800832 0.21893159 -0.300678627
## pct_unemp -0.12862329 0.33009857 -0.030506530
## pct_limited_eng -0.22333346 0.35547555 -0.125707430
## pct_hh_pov -0.30045944 0.24950766 -0.011486308
## dist_m_mine_well cvd_rate_adj res_rate_adj
## mean_pm -0.25751004 0.136954783 0.113962827
## mean_o3 -0.05823036 0.006789463 -0.001793836
## pct_tree_cover 0.47797620 -0.176295758 -0.034307854
## pct_impervious 0.02524829 0.383711656 0.402512635
## mean_aadt_intensity 0.12008641 0.055274160 0.121994187
## dist_m_tri 0.18476022 -0.426652406 -0.376750794
## dist_m_npl 0.03659688 -0.205256148 -0.176011297
## dist_m_waste_site 0.28275484 -0.412987865 -0.424322872
## dist_m_major_emit -0.27042332 0.113286594 -0.012566704
## dist_m_cafo 0.65237893 -0.151743889 -0.061257230
## dist_m_mine_well 1.00000000 -0.410752544 -0.234765650
## cvd_rate_adj -0.41075254 1.000000000 0.864359590
## res_rate_adj -0.23476565 0.864359590 1.000000000
## violent_crime_rate -0.15245003 0.475234675 0.499449246
## property_crime_rate -0.05640681 0.211613232 0.280786581
## pct_less_hs -0.37449548 0.570813439 0.470409304
## pct_unemp -0.27052616 0.468277441 0.441149256
## pct_limited_eng -0.27607853 0.490577454 0.380164971
## pct_hh_pov -0.21843600 0.538480631 0.557944498
## violent_crime_rate property_crime_rate pct_less_hs
## mean_pm 0.05873094 0.088879773 0.14204622
## mean_o3 0.02977242 -0.025156291 0.01633882
## pct_tree_cover -0.20552631 -0.201827442 -0.19081045
## pct_impervious 0.31693831 0.246739067 0.33564418
## mean_aadt_intensity 0.11041575 0.074792588 -0.03883608
## dist_m_tri -0.06251962 -0.058223492 -0.37937630
## dist_m_npl -0.03185242 -0.152322474 -0.15916230
## dist_m_waste_site -0.25287368 -0.268664603 -0.38800832
## dist_m_major_emit -0.01789622 -0.209766780 0.21893159
## dist_m_cafo 0.01675556 -0.064647236 -0.30067863
## dist_m_mine_well -0.15245003 -0.056406808 -0.37449548
## cvd_rate_adj 0.47523468 0.211613232 0.57081344
## res_rate_adj 0.49944925 0.280786581 0.47040930
## violent_crime_rate 1.00000000 0.705151942 0.27967307
## property_crime_rate 0.70515194 1.000000000 -0.00684836
## pct_less_hs 0.27967307 -0.006848360 1.00000000
## pct_unemp 0.35971778 0.007623781 0.62343462
## pct_limited_eng 0.24385889 -0.050184228 0.80738433
## pct_hh_pov 0.42127121 0.163143151 0.72420883
## pct_unemp pct_limited_eng pct_hh_pov
## mean_pm 0.032317153 0.08339659 0.10685820
## mean_o3 0.018969666 0.01164707 -0.01810003
## pct_tree_cover -0.021515177 -0.10974662 0.01115617
## pct_impervious 0.392914001 0.38111402 0.56509450
## mean_aadt_intensity 0.147397153 0.03855846 0.19888999
## dist_m_tri -0.198747643 -0.25490796 -0.31558789
## dist_m_npl 0.132536906 0.01801953 -0.03111065
## dist_m_waste_site -0.128623290 -0.22333346 -0.30045944
## dist_m_major_emit 0.330098571 0.35547555 0.24950766
## dist_m_cafo -0.030506530 -0.12570743 -0.01148631
## dist_m_mine_well -0.270526163 -0.27607853 -0.21843600
## cvd_rate_adj 0.468277441 0.49057745 0.53848063
## res_rate_adj 0.441149256 0.38016497 0.55794450
## violent_crime_rate 0.359717785 0.24385889 0.42127121
## property_crime_rate 0.007623781 -0.05018423 0.16314315
## pct_less_hs 0.623434625 0.80738433 0.72420883
## pct_unemp 1.000000000 0.61293958 0.68048090
## pct_limited_eng 0.612939575 1.00000000 0.74208323
## pct_hh_pov 0.680480902 0.74208323 1.00000000ggplot(pivot_longer(as.data.frame(X.scaled), mean_pm:pct_hh_pov,
names_to = "exp", values_to = "value")) +
geom_histogram(aes(x = value)) +
facet_wrap(~ exp, scales = "free")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
1.2 Covariate data
Covariates were assessed at the individual level. These were selected based on previous HS studies and others in the literature and informed by a DAG.
NOTE: It’ll be interesting to see what comes out of our BEAMERS discussion re: adjusting for gestational age. It’s currently in the analysis
There are four continuous covariates; all of the others have been coded as dummy variables. For the dummy variables, the reference groups are: white_re, ed_grad, norm_bmi
W <- select(hs_data2, latina_re, black_re, other_re,
ed_no_hs, ed_hs, ed_aa, ed_4yr,
low_bmi, ovwt_bmi, obese_bmi,
maternal_age, any_smoker, smokeSH, mean_cpss, mean_epsd,
male, gest_age_w) %>%
as.matrix()
head(W)## latina_re black_re other_re ed_no_hs ed_hs ed_aa ed_4yr low_bmi ovwt_bmi
## [1,] 1 0 0 0 0 1 0 0 0
## [2,] 0 0 1 0 0 1 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 1 0 0 0
## [5,] 0 1 0 0 0 0 1 0 0
## [6,] 1 0 0 0 0 1 0 0 0
## obese_bmi maternal_age any_smoker smokeSH mean_cpss mean_epsd male
## [1,] 0 19 0 1 29 0 0
## [2,] 0 36 0 0 19 2 1
## [3,] 0 34 0 0 19 1 0
## [4,] 0 28 0 0 20 0 0
## [5,] 0 30 0 0 15 0 1
## [6,] 0 22 0 0 17 1 0
## gest_age_w
## [1,] 40.57143
## [2,] 35.85714
## [3,] 40.42857
## [4,] 36.28571
## [5,] 38.42857
## [6,] 40.71429Scaled the non-binary (continuous) covariates
W.s <- apply(W[,c(11, 14, 15, 17)], 2, scale) #' just the continuous ones
W.scaled <- cbind(W[,1:10], W.s[,1],
W[,12:13], W.s[,2:3],
W[,16], W.s[,4])
colnames(W.scaled)## [1] "latina_re" "black_re" "other_re" "ed_no_hs" "ed_hs"
## [6] "ed_aa" "ed_4yr" "low_bmi" "ovwt_bmi" "obese_bmi"
## [11] "" "any_smoker" "smokeSH" "mean_cpss" "mean_epsd"
## [16] "" ""colnames(W.scaled) <- colnames(W)
head(W.scaled)## latina_re black_re other_re ed_no_hs ed_hs ed_aa ed_4yr low_bmi ovwt_bmi
## [1,] 1 0 0 0 0 1 0 0 0
## [2,] 0 0 1 0 0 1 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 1 0 0 0
## [5,] 0 1 0 0 0 0 1 0 0
## [6,] 1 0 0 0 0 1 0 0 0
## obese_bmi maternal_age any_smoker smokeSH mean_cpss mean_epsd male
## [1,] 0 -1.39815187 0 1 3.3147856 -1.2832098 0
## [2,] 0 1.35109608 0 0 0.1179652 -0.6860171 1
## [3,] 0 1.02765515 0 0 0.1179652 -0.9846134 0
## [4,] 0 0.05733234 0 0 0.4376472 -1.2832098 0
## [5,] 0 0.38077328 0 0 -1.1607630 -1.2832098 1
## [6,] 0 -0.91299047 0 0 -0.5213989 -0.9846134 0
## gest_age_w
## [1,] 0.7037686
## [2,] -1.9146645
## [3,] 0.6244221
## [4,] -1.6766251
## [5,] -0.4864283
## [6,] 0.7831150summary(W.scaled)## latina_re black_re other_re ed_no_hs
## Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.0000
## Median :0.0000 Median :0.0000 Median :0.00000 Median :0.0000
## Mean :0.2653 Mean :0.1717 Mean :0.06689 Mean :0.1527
## 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.0000 Max. :1.00000 Max. :1.0000
## ed_hs ed_aa ed_4yr low_bmi
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000
## Median :0.0000 Median :0.0000 Median :0.0000 Median :0.00000
## Mean :0.1851 Mean :0.2319 Mean :0.2185 Mean :0.03344
## 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.00000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.00000
## ovwt_bmi obese_bmi maternal_age any_smoker
## Min. :0.000 Min. :0.0000 Min. :-1.88331 Min. :0.00000
## 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.:-0.91299 1st Qu.:0.00000
## Median :0.000 Median :0.0000 Median : 0.05733 Median :0.00000
## Mean :0.262 Mean :0.1996 Mean : 0.00000 Mean :0.08696
## 3rd Qu.:1.000 3rd Qu.:0.0000 3rd Qu.: 0.70421 3rd Qu.:0.00000
## Max. :1.000 Max. :1.0000 Max. : 2.64486 Max. :1.00000
## smokeSH mean_cpss mean_epsd male
## Min. :0.0000 Min. :-5.9560 Min. :-1.2832 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:-0.5214 1st Qu.:-0.7855 1st Qu.:0.0000
## Median :0.0000 Median : 0.0114 Median :-0.1884 Median :1.0000
## Mean :0.2575 Mean : 0.0000 Mean : 0.0000 Mean :0.5117
## 3rd Qu.:1.0000 3rd Qu.: 0.5442 3rd Qu.: 0.6079 3rd Qu.:1.0000
## Max. :1.0000 Max. : 4.5935 Max. : 6.0324 Max. :1.0000
## gest_age_w
## Min. :-7.7070
## 1st Qu.:-0.3277
## Median : 0.1483
## Mean : 0.0000
## 3rd Qu.: 0.6244
## Max. : 2.9255Variance and histograms for the scaled covariates
var(W.scaled)## latina_re black_re other_re ed_no_hs
## latina_re 0.19514701784 -0.0456034002 -0.0177675585 0.039787884
## black_re -0.04560340022 0.1423669175 -0.0114966555 0.012811803
## other_re -0.01776755853 -0.0114966555 0.0624850693 -0.003531116
## ed_no_hs 0.03978788422 0.0128118032 -0.0035311156 0.129548893
## ed_hs 0.02896808807 0.0150675864 0.0065807155 -0.028296206
## ed_aa 0.01318258282 0.0192967133 0.0023291925 -0.035455487
## ed_4yr -0.03571926262 -0.0163503842 -0.0034713927 -0.033409978
## low_bmi 0.00004479216 0.0020641722 -0.0011235368 -0.003997701
## ovwt_bmi 0.02081218148 -0.0003857103 -0.0052668120 0.004584976
## obese_bmi 0.02065416468 0.0069962872 0.0022620043 0.006318184
## maternal_age -0.10858833047 -0.0895968722 -0.0133074822 -0.137612465
## any_smoker -0.00858889752 0.0185364907 -0.0013586957 0.017954193
## smokeSH -0.00590510033 0.0349789477 0.0084246596 0.030936455
## mean_cpss -0.04668314421 -0.0233123370 0.0255604593 -0.056115354
## mean_epsd 0.04322362524 0.0220578552 0.0202679881 0.067010728
## male -0.00087717989 0.0002202281 0.0003322086 -0.013508570
## gest_age_w -0.02639363168 -0.0368471269 -0.0054750903 -0.007692848
## ed_hs ed_aa ed_4yr low_bmi
## latina_re 0.02896808807 0.013182583 -0.035719263 0.00004479216
## black_re 0.01506758640 0.019296713 -0.016350384 0.00206417224
## other_re 0.00658071548 0.002329193 -0.003471393 -0.00112353679
## ed_no_hs -0.02829620561 -0.035455487 -0.033409978 -0.00399770067
## ed_hs 0.15098194378 -0.042960663 -0.040482163 -0.00173196369
## ed_aa -0.04296066253 0.178312629 -0.050724638 0.00786102484
## ed_4yr -0.04048216276 -0.050724638 0.170951784 -0.00173569637
## low_bmi -0.00173196369 0.007861025 -0.001735696 0.03236233875
## ovwt_bmi -0.00612657270 0.015075052 0.004074843 -0.00877179885
## obese_bmi 0.02441297380 0.009478520 -0.012402453 -0.00668149785
## maternal_age -0.10868310364 -0.040296710 0.109104452 -0.01089529356
## any_smoker 0.00174689441 0.011063665 -0.015673525 0.00155279503
## smokeSH 0.02483352246 0.022806677 -0.036244326 0.00142215122
## mean_cpss -0.03433895561 0.030714793 0.028242932 0.00484169668
## mean_epsd 0.02014619096 0.025424770 -0.046368432 0.00946748435
## male 0.00006345557 0.000630823 0.006367953 -0.00262407430
## gest_age_w -0.01723767893 -0.035168407 0.030147649 -0.00601412877
## ovwt_bmi obese_bmi maternal_age any_smoker smokeSH
## latina_re 0.0208121815 0.020654165 -0.108588330 -0.008588898 -0.005905100
## black_re -0.0003857103 0.006996287 -0.089596872 0.018536491 0.034978948
## other_re -0.0052668120 0.002262004 -0.013307482 -0.001358696 0.008424660
## ed_no_hs 0.0045849757 0.006318184 -0.137612465 0.017954193 0.030936455
## ed_hs -0.0061265727 0.024412974 -0.108683104 0.001746894 0.024833522
## ed_aa 0.0150750518 0.009478520 -0.040296710 0.011063665 0.022806677
## ed_4yr 0.0040748427 -0.012402453 0.109104452 -0.015673525 -0.036244326
## low_bmi -0.0087717989 -0.006681498 -0.010895294 0.001552795 0.001422151
## ovwt_bmi 0.1935643614 -0.052338400 0.008900228 -0.006065606 -0.010623208
## obese_bmi -0.0523383998 0.159910515 0.002790074 0.002717391 0.011052467
## maternal_age 0.0089002276 0.002790074 1.000000000 -0.046629611 -0.155964054
## any_smoker -0.0060656056 0.002717391 -0.046629611 0.079483696 0.049010093
## smokeSH -0.0106232083 0.011052467 -0.155964054 0.049010093 0.191419314
## mean_cpss -0.0082476896 -0.008722611 0.100637638 0.017642908 0.031721118
## mean_epsd -0.0019127169 0.028133038 -0.160410684 0.042144665 0.108180210
## male 0.0008361204 -0.001780489 0.023413804 0.002329193 0.002004449
## gest_age_w -0.0148466585 -0.021461430 0.091663607 -0.014981418 -0.050311537
## mean_cpss mean_epsd male gest_age_w
## latina_re -0.046683144 0.043223625 -0.00087717989 -0.026393632
## black_re -0.023312337 0.022057855 0.00022022814 -0.036847127
## other_re 0.025560459 0.020267988 0.00033220855 -0.005475090
## ed_no_hs -0.056115354 0.067010728 -0.01350857023 -0.007692848
## ed_hs -0.034338956 0.020146191 0.00006345557 -0.017237679
## ed_aa 0.030714793 0.025424770 0.00063082298 -0.035168407
## ed_4yr 0.028242932 -0.046368432 0.00636795270 0.030147649
## low_bmi 0.004841697 0.009467484 -0.00262407430 -0.006014129
## ovwt_bmi -0.008247690 -0.001912717 0.00083612040 -0.014846658
## obese_bmi -0.008722611 0.028133038 -0.00178048853 -0.021461430
## maternal_age 0.100637638 -0.160410684 0.02341380415 0.091663607
## any_smoker 0.017642908 0.042144665 0.00232919255 -0.014981418
## smokeSH 0.031721118 0.108180210 0.00200444935 -0.050311537
## mean_cpss 1.000000000 0.455187203 -0.00331530432 -0.037142336
## mean_epsd 0.455187203 1.000000000 0.00154181454 -0.137187808
## male -0.003315304 0.001541815 0.25014184185 -0.007427180
## gest_age_w -0.037142336 -0.137187808 -0.00742717951 1.000000000ggplot(pivot_longer(as.data.frame(W.scaled), latina_re:gest_age_w,
names_to = "exp", values_to = "value")) +
geom_histogram(aes(x = value)) +
facet_wrap(~ exp, scales = "free")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
1.3 Response data: birth weight
Y <- select(hs_data2, birth_weight) %>%
as.matrix()
head(Y)## birth_weight
## [1,] 2860
## [2,] 2755
## [3,] 3505
## [4,] 2695
## [5,] 3355
## [6,] 3810Distribution of birth weight and scaled birth weight
hist(Y, breaks = 20)
hist(scale(Y), breaks = 20)
1.4 Scatterplots of exposures and outcome (birth weight)
Both birth weight (Y) and the exposures are scaled here
NOTE: Don’t use these plots as a way to estimate how many predictors might make the cut. This should be done a priori
df <- as.data.frame(cbind(scale(Y), X.scaled))
# par(mfrow=c(5,4))
sapply(2:length(df), function(x){
lm.x <- lm(birth_weight ~ df[,x], data = df)
plot(df[,c(x, 1)],
xlab = paste0(colnames(df)[x], " beta: ",
round(summary(lm.x)$coef[2,1],4),
"; p = ",
round(summary(lm.x)$coef[2,4],4)))
abline(lm.x)
})
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
##
## [[7]]
## NULL
##
## [[8]]
## NULL
##
## [[9]]
## NULL
##
## [[10]]
## NULL
##
## [[11]]
## NULL
##
## [[12]]
## NULL
##
## [[13]]
## NULL
##
## [[14]]
## NULL
##
## [[15]]
## NULL
##
## [[16]]
## NULL
##
## [[17]]
## NULL
##
## [[18]]
## NULL
##
## [[19]]
## NULL2 Is gestational age a potenital mediator?
I.e., is there a relationship between our exposures and gestational age?
The DAG might look something like this:
exposures —> gestational age —> birth weight
2.1 Scatter plots for exposures and gestational age
Both gestational age and the exposures are scaled here. Gestational age measured in weeks from estimated date of conception to delivery
Since there were some (small) relationships between exposures and gestational age (based on simple linear regression models– namely the ozone and SES indicators), I’m going to omit this covariate for now.
df2 <- as.data.frame(cbind(W.scaled[,"gest_age_w"], X.scaled))
colnames(df2)[1] <- "gest_age_w"
# par(mfrow=c(5,4))
sapply(2:length(df2), function(x){
lm.x <- lm(gest_age_w ~ df2[,x], data = df2)
plot(df2[,c(x, 1)],
xlab = paste0(colnames(df2)[x], " beta: ",
round(summary(lm.x)$coef[2,1],4),
"; p = ",
round(summary(lm.x)$coef[2,4],4)))
abline(lm.x)
})
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
##
## [[7]]
## NULL
##
## [[8]]
## NULL
##
## [[9]]
## NULL
##
## [[10]]
## NULL
##
## [[11]]
## NULL
##
## [[12]]
## NULL
##
## [[13]]
## NULL
##
## [[14]]
## NULL
##
## [[15]]
## NULL
##
## [[16]]
## NULL
##
## [[17]]
## NULL
##
## [[18]]
## NULL
##
## [[19]]
## NULLDropping gest_age_w from the covariates
W.scaled2 <- W.scaled[,-c(ncol(W.scaled))]
colnames(W.scaled2)## [1] "latina_re" "black_re" "other_re" "ed_no_hs" "ed_hs"
## [6] "ed_aa" "ed_4yr" "low_bmi" "ovwt_bmi" "obese_bmi"
## [11] "maternal_age" "any_smoker" "smokeSH" "mean_cpss" "mean_epsd"
## [16] "male"3 RIDGE regression
To see if there might be something going on, Lauren suggested a ridge regression with a small penalty.
set.seed(123)
library(glmnet)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 3.0-2lambda_seq <- 10^seq(2, -2, by = -.1)
#' Best lambda from CV
ridge_cv <- cv.glmnet(X.scaled, scale(Y), alpha = 0, lambda = lambda_seq)
plot(ridge_cv)
best_lambda <- ridge_cv$lambda.min
best_lambda## [1] 1.258925#' Fit the model using the best_lambda
bw_ridge <- glmnet(X.scaled, scale(Y), alpha = 0, lambda = best_lambda)
summary(bw_ridge)## Length Class Mode
## a0 1 -none- numeric
## beta 19 dgCMatrix S4
## df 1 -none- numeric
## dim 2 -none- numeric
## lambda 1 -none- numeric
## dev.ratio 1 -none- numeric
## nulldev 1 -none- numeric
## npasses 1 -none- numeric
## jerr 1 -none- numeric
## offset 1 -none- logical
## call 5 -none- call
## nobs 1 -none- numericRidge regression coefficients
coef(bw_ridge)## 20 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) -0.000000000000000214871
## mean_pm 0.009001292934592152947
## mean_o3 -0.048621510737574098748
## pct_tree_cover 0.001342271308724626074
## pct_impervious -0.012521694754762259516
## mean_aadt_intensity -0.005312589097233479454
## dist_m_tri -0.001659989297196438852
## dist_m_npl 0.001245550042200224651
## dist_m_waste_site 0.008057207316642686287
## dist_m_major_emit -0.002818043274053380360
## dist_m_cafo -0.003389166474688925200
## dist_m_mine_well -0.009950869706335926587
## cvd_rate_adj -0.016027916742713192721
## res_rate_adj -0.009686648757473777238
## violent_crime_rate -0.005457924031598470650
## property_crime_rate -0.000364961614232602028
## pct_less_hs -0.021599485837649336911
## pct_unemp -0.042979714186389718356
## pct_limited_eng -0.010315534078569870563
## pct_hh_pov -0.010927666957598947128Ridge regression predictions
ridge_pred <- predict(bw_ridge, newx = X.scaled)
plot(scale(Y), ridge_pred)
actual <- scale(Y)
preds <- ridge_pred
rsq <- 1 - (sum((preds - actual) ^ 2))/(sum((actual - mean(actual)) ^ 2))The R2 value for this model is 0.03. Based on these results, it doesn’t look like there’s much here.
4 Nonparametric Bayesian Shrinkage (NPB): Birth weight
Still, we wanted to try to fit the NPB model with these data.
4.1 Finding the NPB priors
Start with Lauren’s from the example in the vignette
In an email from April 29, Lauren provided me with some additional guidance on finding the NPB priors:
- Keep alpha.pi and beta.pi set to 1, and then let a.phi1 take values 1, 10, and 100 and see how the results change.
- Keep a.phi1 set to 1 (or 10 or 100), and mess with alpha.pi and beta.pi.
- Run the following code: alpha.pi=1 beta.pi=1 plot(density(rbeta(10000, alpha.pi, beta.pi))) and then change alpha.pi and beta.pi and see how it changes the prior distribution.
- This is the distribution of the probability of a main effect regression coefficient being 0 (aka exclusion probability). We don’t want this to be too informative (you don’t want high mass around just a few values). Also alpha.pi and beta.pi don’t have to be the same value. You might try alpha.pi = 1 and beta.pi = 2 to get a slightly lower prior probability of exclusion. Try not to change all three (alpha.pi, beta.pi, and a.phi1) at once.
- When playing with the priors, set “interact=FALSE” and just fit the modelwith the main effects. Most of the interactions were null anyway so it shouldn’t change the results too much and it will make the code run a lot faster. Then when you find a set of priors you like, you can add in “interact=TRUE” and “XWinteract=TRUE.”
Some additional feedback from Lauren during our 6/10 meeting:
- The confidence intervals were really wide and heavily skewed. I’m going to try adjusting the sig2inv.mu1 parameter after the a.phi1 parameter to see if this helps
- the rbeta distributions is interpreted as the exclusion probability, so I should try to aim to have most of the mass of that distribution in the middle, since we hypothesize that maybe 40-60% of the predictors will be important. The way to do this is to set alpha.pi and beta.pi to the same value
- I should set the burn in number to be about half the iterations
4.1.1 Vignette Priors
set.seed(123)
priors.npb.1 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1)
fit.npb.1 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.1, interact = F)
npb.sum.1 <- summary(fit.npb.1)
npb.sum.1$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.42293623 3.1955064 0.000000 0 0.022
## [2,] -7.77734425 17.0869661 -61.452320 0 0.208
## [3,] -0.34921645 2.9046941 0.000000 0 0.020
## [4,] -0.72527602 4.5329148 -9.113312 0 0.032
## [5,] -0.10902198 1.6231693 0.000000 0 0.008
## [6,] -0.22439232 3.4307474 0.000000 0 0.022
## [7,] 0.13813603 3.2405735 0.000000 0 0.014
## [8,] 0.50835031 5.5683575 0.000000 0 0.018
## [9,] -0.01225732 0.1949014 0.000000 0 0.004
## [10,] -0.08255604 1.1888837 0.000000 0 0.006
## [11,] -1.04104994 5.8717128 -21.539810 0 0.050
## [12,] -0.44293941 3.6077235 0.000000 0 0.020
## [13,] -1.34313104 6.6180314 -26.725959 0 0.050
## [14,] -0.21174177 1.9565687 0.000000 0 0.016
## [15,] -0.67619819 4.4923381 -7.070728 0 0.034
## [16,] -0.69671736 4.5842905 -3.417542 0 0.026
## [17,] -3.11238035 10.7411492 -42.631782 0 0.096
## [18,] -0.58256651 4.1194858 -3.681873 0 0.032
## [19,] -0.41969750 3.4549837 0.000000 0 0.020plot(fit.npb.1$beta[,1], type = "l")
plot(fit.npb.1$beta[,2], type = "l")
plot(fit.npb.1$beta[,13], type = "l")
4.1.2 Adjust sig2inv.mu1
The default value for this parameter is 1. Adjusting this parameter may help narrow these confidence intervals a bit
4.1.2.1 Try sig2inv.mu1 = 10
That seemed to make the confidence intervals wider
priors.npb.2 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
sig2inv.mu1 = 10)
fit.npb.2 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.2, interact = F)
npb.sum.2 <- summary(fit.npb.2)
npb.sum.2$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.57320485 3.542448 -10.457164 0.000000 0.060
## [2,] -9.76615860 18.286535 -58.159264 0.000000 0.294
## [3,] -0.48674498 2.686788 -8.315099 0.000000 0.056
## [4,] -1.17031436 6.050804 -16.685042 0.000000 0.084
## [5,] -0.18426760 3.095979 -4.889139 0.000000 0.042
## [6,] -0.17137354 3.496180 -7.318933 0.000000 0.058
## [7,] 0.12313460 3.434666 -1.841864 0.000000 0.054
## [8,] 0.68066167 6.178553 0.000000 7.296187 0.050
## [9,] 0.07551352 3.152292 0.000000 0.000000 0.038
## [10,] -0.24556275 1.952749 -6.237070 0.000000 0.048
## [11,] -1.22843970 5.154892 -17.980677 0.000000 0.096
## [12,] -0.62684355 3.672637 -9.382015 0.000000 0.058
## [13,] -2.42125423 8.397953 -34.673034 0.000000 0.122
## [14,] -0.43994939 3.003203 -8.243475 0.000000 0.052
## [15,] -1.12213230 4.786993 -16.736005 0.000000 0.078
## [16,] -1.20895702 5.295603 -18.793484 0.000000 0.086
## [17,] -4.07327965 11.899635 -49.929697 0.000000 0.158
## [18,] -0.47703166 2.789470 -8.900533 0.000000 0.064
## [19,] -0.81382810 4.258625 -11.526663 0.000000 0.070plot(fit.npb.2$beta[,1], type = "l")
plot(fit.npb.2$beta[,2], type = "l")
plot(fit.npb.2$beta[,13], type = "l")
4.1.2.2 Try sig2inv.mu1 = 100
That also seemed to make the confidence intervals wider
priors.npb.3 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
sig2inv.mu1 = 100)
fit.npb.3 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.3, interact = F)
npb.sum.3 <- summary(fit.npb.3)
npb.sum.3$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.095131290 3.140375 0.000000 0.0000000 0.020
## [2,] -11.681887793 20.022155 -59.403940 0.0000000 0.292
## [3,] 0.016920655 3.934123 0.000000 0.0000000 0.028
## [4,] -0.914938821 5.470093 -18.182994 0.0000000 0.040
## [5,] -0.248071932 2.846108 0.000000 0.0000000 0.034
## [6,] -0.059851371 2.346589 0.000000 0.0000000 0.030
## [7,] 0.517100918 4.544226 0.000000 0.2186056 0.032
## [8,] 1.371237081 8.156704 0.000000 17.0215314 0.050
## [9,] 0.252087787 2.563337 0.000000 0.0000000 0.026
## [10,] -0.003136574 2.221280 0.000000 0.0000000 0.030
## [11,] -0.871617648 5.610285 -18.258161 0.0000000 0.040
## [12,] -0.605872070 4.903910 -12.928430 0.0000000 0.044
## [13,] -1.138512749 6.681420 -23.502834 0.0000000 0.058
## [14,] -0.053623560 2.120945 0.000000 0.0000000 0.024
## [15,] -0.430496025 4.375130 -2.256698 0.0000000 0.042
## [16,] -0.556105327 7.871809 -22.192950 0.0000000 0.058
## [17,] -4.035417542 12.806278 -49.923011 0.0000000 0.120
## [18,] -0.334067046 3.709821 0.000000 0.0000000 0.036
## [19,] -1.149915311 7.331426 -10.142257 0.0000000 0.030plot(fit.npb.3$beta[,1], type = "l")
plot(fit.npb.3$beta[,2], type = "l")
plot(fit.npb.3$beta[,13], type = "l")
4.1.2.3 Try sig2inv.mu1 = 0.1?
This might be too narrow?
priors.npb.4 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
sig2inv.mu1 = 0.1)
fit.npb.4 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.4, interact = F)
npb.sum.4 <- summary(fit.npb.4)
npb.sum.4$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.105976545 1.92589731 0.00000 0 0.008
## [2,] -3.708253022 12.88063652 -51.92299 0 0.086
## [3,] -0.004328270 0.09383458 0.00000 0 0.004
## [4,] -0.133231949 2.20638272 0.00000 0 0.006
## [5,] -0.010487186 0.16648611 0.00000 0 0.006
## [6,] -0.003353958 0.07499678 0.00000 0 0.002
## [7,] 0.000000000 0.00000000 0.00000 0 0.000
## [8,] 0.226086189 3.07843115 0.00000 0 0.010
## [9,] -0.015321283 1.31966148 0.00000 0 0.006
## [10,] -0.127552312 1.86015024 0.00000 0 0.008
## [11,] -0.385794753 4.12927001 0.00000 0 0.010
## [12,] -0.400850294 4.05400654 0.00000 0 0.014
## [13,] -0.786528947 5.60881484 0.00000 0 0.024
## [14,] -0.112417701 1.86716884 0.00000 0 0.004
## [15,] -0.266540085 3.04248660 0.00000 0 0.012
## [16,] -0.063191855 1.00336679 0.00000 0 0.004
## [17,] -1.612479624 8.44032992 -36.35425 0 0.040
## [18,] -0.098734619 1.42912623 0.00000 0 0.008
## [19,] -0.238959635 3.10000213 0.00000 0 0.008plot(fit.npb.4$beta[,1], type = "l")
plot(fit.npb.4$beta[,2], type = "l")
plot(fit.npb.4$beta[,13], type = "l")
4.1.2.4 Try sig2inv.mu1 = 0.5?
Also too narrow– going back to 1
priors.npb.4 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
sig2inv.mu1 = 0.1)
fit.npb.4 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.4, interact = F)
npb.sum.4 <- summary(fit.npb.4)
npb.sum.4$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.16653792 2.1552452 0.000000 0 0.016
## [2,] -3.07753903 11.8843477 -50.000667 0 0.082
## [3,] -0.05255064 0.9993058 0.000000 0 0.006
## [4,] -0.28271383 2.5061098 0.000000 0 0.022
## [5,] -0.06295496 1.1247119 0.000000 0 0.006
## [6,] 0.02999439 0.7687866 0.000000 0 0.006
## [7,] 0.21649613 3.4418344 0.000000 0 0.008
## [8,] 0.50724409 4.6522999 0.000000 0 0.014
## [9,] -0.03404257 0.5868220 0.000000 0 0.004
## [10,] -0.12325749 1.8669486 0.000000 0 0.014
## [11,] -0.20457023 2.3663177 0.000000 0 0.010
## [12,] -0.29222742 3.2766185 0.000000 0 0.016
## [13,] -0.41257543 3.4489993 0.000000 0 0.020
## [14,] -0.13817717 1.3110905 0.000000 0 0.012
## [15,] -0.21872498 3.3383078 0.000000 0 0.020
## [16,] -0.24942108 2.5110565 0.000000 0 0.014
## [17,] -0.64575618 4.5204832 -4.510545 0 0.026
## [18,] -0.25473085 2.7429636 0.000000 0 0.016
## [19,] -0.40524422 2.9289043 0.000000 0 0.024plot(fit.npb.4$beta[,1], type = "l")
plot(fit.npb.4$beta[,2], type = "l")
plot(fit.npb.4$beta[,13], type = "l")
4.1.3 Next, try adjusting a.phi1
This parameter is the shape parameter for gamma prior on phi sub1 inv2. It’s default value is 1
4.1.3.1 Try making a.phi1 = 10
PIPs again increase slightly
priors.npb.5 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10)
fit.npb.5 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.5, interact = F)
npb.sum.5 <- summary(fit.npb.5)
npb.sum.5$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.581552391 3.821619 -6.9203902 0.00000 0.040
## [2,] -11.473528795 20.634360 -61.5619258 0.00000 0.282
## [3,] -0.280908392 2.581600 0.0000000 0.00000 0.026
## [4,] -0.685781042 4.763465 -11.8478088 0.00000 0.048
## [5,] 0.007128347 1.470303 0.0000000 0.00000 0.020
## [6,] -0.063487381 2.012826 0.0000000 0.00000 0.026
## [7,] 0.346319628 4.162866 0.0000000 0.00000 0.030
## [8,] 1.383203139 7.315417 0.0000000 28.75221 0.050
## [9,] 0.046533892 1.936927 0.0000000 0.00000 0.032
## [10,] -0.223120797 2.664584 0.0000000 0.00000 0.024
## [11,] -1.097947666 5.698955 -18.5451922 0.00000 0.056
## [12,] -0.584701314 5.176760 -5.6821139 0.00000 0.036
## [13,] -1.339277024 6.544944 -26.5626778 0.00000 0.058
## [14,] -0.093226827 2.408309 0.0000000 0.00000 0.030
## [15,] -0.525808339 3.716516 -5.3152818 0.00000 0.032
## [16,] -0.851948191 5.505950 -13.1334253 0.00000 0.040
## [17,] -4.359826211 13.001580 -50.7884608 0.00000 0.136
## [18,] -0.251327971 2.127834 -0.5059922 0.00000 0.032
## [19,] -0.965748807 6.050453 -14.9412688 0.00000 0.052plot(fit.npb.5$beta[,1], type = "l")
plot(fit.npb.5$beta[,2], type = "l")
plot(fit.npb.5$beta[,13], type = "l")
4.1.3.2 Try making a.phi1 = 10 and sig2inv.mu1 = 10
Definitely worse in terms of the confidence intervals
priors.npb.6 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 10)
fit.npb.6 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.6, interact = F)
npb.sum.6 <- summary(fit.npb.6)
npb.sum.6$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.4928015 3.286066 -4.847052 0.00000 0.034
## [2,] -17.4623243 22.302478 -64.918518 0.00000 0.452
## [3,] -0.2814630 2.585887 -2.677709 0.00000 0.034
## [4,] -2.0754808 8.100446 -33.727567 0.00000 0.088
## [5,] 0.1666961 3.909889 0.000000 0.00000 0.032
## [6,] -0.0959787 3.649590 0.000000 0.00000 0.024
## [7,] 0.2210554 3.615359 0.000000 0.00000 0.040
## [8,] 2.2210034 10.266335 0.000000 40.50969 0.072
## [9,] 0.1968973 2.996566 0.000000 0.00000 0.030
## [10,] -0.1913085 2.891537 -4.718272 0.00000 0.038
## [11,] -1.8156249 7.942042 -29.464176 0.00000 0.094
## [12,] -0.7893540 4.525473 -12.672066 0.00000 0.052
## [13,] -2.9937749 10.303755 -43.520723 0.00000 0.110
## [14,] -0.2091624 3.483150 -8.291385 0.00000 0.050
## [15,] -1.1753419 5.587383 -17.375845 0.00000 0.068
## [16,] -1.5192580 7.132802 -26.391704 0.00000 0.070
## [17,] -7.1958362 15.741948 -53.830439 0.00000 0.218
## [18,] -0.3065605 4.195224 -8.740062 0.00000 0.046
## [19,] -1.0956311 6.006086 -19.708031 0.00000 0.052plot(fit.npb.6$beta[,1], type = "l")
plot(fit.npb.6$beta[,2], type = "l")
plot(fit.npb.6$beta[,13], type = "l")
4.1.3.3 Try making a.phi1 = 10 and sig2inv.mu1 = 100
Not an improvement
priors.npb.7 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 100)
fit.npb.7 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.7, interact = F)
npb.sum.7 <- summary(fit.npb.7)
npb.sum.7$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.22104871 2.757115 0.000000 0.00000 0.038
## [2,] -20.41085141 23.151185 -66.092224 0.00000 0.524
## [3,] -0.24631131 2.661488 0.000000 0.00000 0.034
## [4,] -2.09887813 8.881620 -35.022901 0.00000 0.072
## [5,] 0.11682112 4.046631 0.000000 0.00000 0.032
## [6,] 0.20596362 2.984464 0.000000 0.00000 0.028
## [7,] 0.38192374 3.429390 0.000000 0.00000 0.038
## [8,] 2.55506778 10.355460 0.000000 44.07723 0.082
## [9,] 0.26757819 3.207645 0.000000 0.00000 0.030
## [10,] -0.21437574 3.450039 0.000000 0.00000 0.036
## [11,] -1.88737645 7.711836 -30.378845 0.00000 0.074
## [12,] -1.04124580 7.530258 -21.971510 0.00000 0.084
## [13,] -3.43025706 11.053184 -43.266666 0.00000 0.122
## [14,] 0.03994552 4.879433 -5.159376 0.00000 0.046
## [15,] -1.57907973 7.007126 -23.564577 0.00000 0.084
## [16,] -1.07553220 6.226368 -21.246083 0.00000 0.064
## [17,] -7.66359663 16.478164 -57.107719 0.00000 0.218
## [18,] -0.62459729 4.766040 -14.155347 0.00000 0.050
## [19,] -0.94616666 5.567583 -17.775554 0.00000 0.058plot(fit.npb.7$beta[,1], type = "l")
plot(fit.npb.7$beta[,2], type = "l")
plot(fit.npb.7$beta[,13], type = "l")
4.1.3.4 Try making a.phi1 = 10 and sig2inv.mu1 = 0.1
Nope
priors.npb.8 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 0.1)
fit.npb.8 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.8, interact = F)
npb.sum.8 <- summary(fit.npb.8)
npb.sum.8$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.15754855 2.1361739 0.000000 0 0.016
## [2,] -6.30991048 16.1308085 -56.736167 0 0.160
## [3,] -0.11841473 1.8581095 0.000000 0 0.018
## [4,] -0.34787660 3.0322241 0.000000 0 0.020
## [5,] -0.04711339 0.7135716 0.000000 0 0.008
## [6,] 0.24920300 3.4996638 0.000000 0 0.024
## [7,] 0.01351782 0.7029928 0.000000 0 0.010
## [8,] 0.39877502 3.9120930 0.000000 0 0.028
## [9,] 0.11288992 2.0333384 0.000000 0 0.020
## [10,] -0.09659308 1.3557107 0.000000 0 0.014
## [11,] -0.55699468 4.2665677 0.000000 0 0.026
## [12,] -0.74334034 5.3325922 -9.749636 0 0.044
## [13,] -1.22844989 6.3929916 -23.584655 0 0.044
## [14,] -0.01716576 1.0229667 0.000000 0 0.014
## [15,] -0.35293666 3.0341882 0.000000 0 0.018
## [16,] -0.50600924 3.9284400 0.000000 0 0.030
## [17,] -1.72631176 8.0307307 -34.544925 0 0.058
## [18,] -0.20194757 2.5711228 0.000000 0 0.020
## [19,] -0.27606966 2.9828552 0.000000 0 0.016plot(fit.npb.8$beta[,1], type = "l")
plot(fit.npb.8$beta[,2], type = "l")
plot(fit.npb.8$beta[,13], type = "l")
4.1.3.5 Try making a.phi1 = 100
priors.npb.9 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 1)
fit.npb.9 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.9, interact = F)
npb.sum.9 <- summary(fit.npb.9)
npb.sum.9$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.49081449 4.629401 -11.061650 0.0000000 0.072
## [2,] -14.61898497 21.824480 -64.166897 0.0000000 0.384
## [3,] -0.42735467 3.520918 -9.044472 0.0000000 0.080
## [4,] -1.32439326 7.050419 -24.968932 0.0000000 0.090
## [5,] -0.16896624 3.790224 -5.342744 0.0000000 0.064
## [6,] 0.08680043 3.475061 -1.363339 0.7361937 0.054
## [7,] 0.21923467 4.910740 -5.391763 7.9913316 0.074
## [8,] 3.71035015 12.113573 0.000000 47.9510551 0.142
## [9,] 0.80244892 5.594285 0.000000 15.8078721 0.064
## [10,] -0.01292279 3.611167 -4.312884 0.9968086 0.064
## [11,] -1.62260660 7.209524 -26.960978 0.0000000 0.094
## [12,] -0.50405112 4.702056 -10.830387 0.2037986 0.080
## [13,] -2.37992543 9.107271 -34.744576 0.0000000 0.108
## [14,] 0.14853855 4.109383 -2.460935 0.2037986 0.062
## [15,] -0.61961911 3.762306 -8.832907 0.0000000 0.062
## [16,] -0.96949695 6.278286 -20.260757 0.0000000 0.092
## [17,] -5.73190550 14.774390 -52.281428 0.0000000 0.188
## [18,] -0.82480285 4.947323 -17.562099 0.0000000 0.092
## [19,] -1.39328755 7.002128 -22.543324 0.0000000 0.092plot(fit.npb.9$beta[,1], type = "l")
plot(fit.npb.9$beta[,2], type = "l")
plot(fit.npb.9$beta[,13], type = "l")
4.1.3.6 Try making a.phi1 = 100 and sig2inv.mu1 = 10
priors.npb.10 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 10)
fit.npb.10 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.10, interact = F)
npb.sum.10 <- summary(fit.npb.10)
npb.sum.10$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.8955066 4.608069 -14.166614 2.0820725 0.156
## [2,] -22.3467639 22.833264 -67.216503 0.0000000 0.632
## [3,] -1.0899915 6.266198 -17.583821 4.5550541 0.162
## [4,] -2.3737495 8.805426 -33.820739 4.7554120 0.202
## [5,] -0.3559511 4.337999 -11.343352 6.0657719 0.142
## [6,] 0.1227590 5.098899 -11.856142 13.5502226 0.150
## [7,] 1.0923054 6.925488 -7.915575 22.8299239 0.154
## [8,] 4.2048826 13.169188 -6.612679 45.9947819 0.230
## [9,] 0.3450384 4.113284 -7.279409 11.8604506 0.138
## [10,] -0.3254172 4.709016 -10.877562 7.9023656 0.140
## [11,] -3.5999901 9.997178 -34.815423 0.6276695 0.224
## [12,] -1.5518494 8.235853 -23.908544 6.9765096 0.192
## [13,] -3.6389730 10.311734 -37.952368 2.9627703 0.252
## [14,] -0.2175783 5.652138 -13.317473 11.8566847 0.158
## [15,] -1.6873238 6.697213 -22.019093 1.7058149 0.162
## [16,] -2.5428487 9.312482 -32.166143 6.5707112 0.210
## [17,] -8.0393071 16.401686 -53.826942 0.0000000 0.310
## [18,] -1.1660232 6.244023 -20.317334 6.9765096 0.176
## [19,] -1.7402131 7.377813 -26.084992 3.0057293 0.182plot(fit.npb.10$beta[,1], type = "l")
plot(fit.npb.10$beta[,2], type = "l")
plot(fit.npb.10$beta[,13], type = "l")
4.1.3.7 Try making a.phi1 = 100 and sig2inv.mu1 = 100
Still not great.
priors.npb.11 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 100)
fit.npb.11 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.11, interact = F)
npb.sum.11 <- summary(fit.npb.11)
npb.sum.11$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.7950218 5.244945 -16.970623 4.537841 0.126
## [2,] -27.2058531 21.929493 -68.156952 0.000000 0.728
## [3,] -0.4288572 4.965463 -13.452307 3.912050 0.116
## [4,] -3.2572407 9.847259 -36.272125 1.027606 0.196
## [5,] 0.2192220 4.831992 -5.568812 12.263959 0.132
## [6,] 0.2665736 6.216597 -9.888779 14.244300 0.128
## [7,] 1.0229973 6.669657 -4.270251 21.779034 0.148
## [8,] 3.9530248 11.917094 -2.002081 44.514415 0.198
## [9,] 1.2600757 7.541149 -8.009803 27.325068 0.172
## [10,] -0.1415800 6.526332 -12.783186 8.690473 0.134
## [11,] -3.3426308 10.322144 -38.047100 0.000000 0.172
## [12,] -1.2593536 8.369107 -24.126005 9.884446 0.164
## [13,] -3.5622279 11.517625 -40.880095 5.691475 0.212
## [14,] 0.5091027 6.421518 -13.023474 23.268069 0.156
## [15,] -2.4704408 8.688523 -31.742324 4.421180 0.198
## [16,] -2.6022436 10.254831 -36.585588 5.991238 0.182
## [17,] -10.5991456 17.509123 -51.783170 0.000000 0.362
## [18,] -1.4959087 7.371900 -28.144047 4.255691 0.142
## [19,] -2.1099624 8.284785 -28.971796 3.879175 0.160plot(fit.npb.11$beta[,1], type = "l")
plot(fit.npb.11$beta[,2], type = "l")
plot(fit.npb.11$beta[,13], type = "l")
4.1.4 Adjust alpha.pi and beta.pi
For now, leave a.phi1 and sig2inv.mu1 alone for now.
alpha.pi and beta.pi are responisble for the exclusion probability distribution. If we thing we want ~50% of our covariates, we need the mass of this distribution to be somewhere between 0.4 and 0.6. To do this, we set alpha.pi and beta.pi to the same value
4.1.4.1 Try making alpha.pi and beta.pi 2
plot(density(rbeta(10000, 2, 2)))
priors.npb.12 <- list(alpha.pi = 2, beta.pi = 2, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 1)
fit.npb.12 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.12, interact = F)
npb.sum.12 <- summary(fit.npb.12)
npb.sum.12$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.548121233 4.0985444 -8.338155 0.00000 0.040
## [2,] -11.687879813 19.4861619 -59.865735 0.00000 0.330
## [3,] -0.377632766 3.5324695 -1.735183 0.00000 0.034
## [4,] -1.194321283 5.7745675 -22.468250 0.00000 0.060
## [5,] -0.054825515 0.8274353 0.000000 0.00000 0.010
## [6,] 0.062375455 2.4983039 0.000000 0.00000 0.028
## [7,] -0.036286670 2.1712654 0.000000 0.00000 0.024
## [8,] 1.280627308 7.7689013 0.000000 31.03455 0.048
## [9,] 0.008804336 0.6536777 0.000000 0.00000 0.008
## [10,] -0.196413518 1.6544415 0.000000 0.00000 0.020
## [11,] -1.760859372 7.2421988 -28.993435 0.00000 0.072
## [12,] -0.836066000 4.7927503 -16.370721 0.00000 0.056
## [13,] -1.918927858 7.8922310 -30.315590 0.00000 0.084
## [14,] -0.316961255 2.2995747 0.000000 0.00000 0.024
## [15,] -0.842456591 5.1178116 -13.396083 0.00000 0.038
## [16,] -1.653424484 6.8653248 -27.957693 0.00000 0.070
## [17,] -5.446856463 13.2720490 -46.601885 0.00000 0.180
## [18,] -0.512403416 3.4378903 -6.786906 0.00000 0.038
## [19,] -0.772038898 4.4605108 -13.288182 0.00000 0.044plot(fit.npb.12$beta[,1], type = "l")
plot(fit.npb.12$beta[,2], type = "l")
plot(fit.npb.12$beta[,13], type = "l")
4.1.4.2 Try making alpha.pi and beta.pi 5
plot(density(rbeta(10000, 5, 5)))
priors.npb.13 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 1)
fit.npb.13 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.13, interact = F)
npb.sum.13 <- summary(fit.npb.13)
npb.sum.13$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.5175818 6.596365 -24.41351 0.0000000 0.224
## [2,] -17.5507282 19.138119 -62.21890 0.0000000 0.658
## [3,] -1.3631356 4.765574 -16.61580 0.0000000 0.172
## [4,] -2.8709572 7.662579 -26.46499 0.0000000 0.234
## [5,] -1.0397884 4.671120 -15.15083 0.6703999 0.138
## [6,] -0.9327580 5.341849 -16.01600 3.3862704 0.168
## [7,] -0.4536953 4.155734 -11.29585 1.9789806 0.112
## [8,] 0.7424446 8.219155 -10.60626 28.5344189 0.136
## [9,] -0.3941859 4.179095 -11.10859 4.1061592 0.132
## [10,] -0.7716985 4.078803 -12.39498 2.0461425 0.128
## [11,] -4.6362502 10.539032 -39.73617 0.4258184 0.284
## [12,] -2.4914016 7.135969 -27.08012 0.0000000 0.224
## [13,] -4.9526503 10.546835 -38.60816 0.0000000 0.302
## [14,] -0.9026749 4.814544 -14.34211 1.3585607 0.150
## [15,] -3.1748864 8.047856 -31.50061 0.0000000 0.238
## [16,] -2.8539683 7.600540 -25.40936 1.1146465 0.252
## [17,] -9.3501374 15.413350 -51.29416 0.0000000 0.420
## [18,] -1.8299513 5.819704 -18.76603 0.0000000 0.188
## [19,] -2.3426401 6.320223 -24.31243 0.0000000 0.218plot(fit.npb.13$beta[,1], type = "l")
plot(fit.npb.13$beta[,2], type = "l")
plot(fit.npb.13$beta[,13], type = "l")
4.1.4.3 Try making alpha.pi and beta.pi 8
plot(density(rbeta(10000, 8, 8)))
priors.npb.14 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 1)
fit.npb.14 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.14, interact = F)
npb.sum.14 <- summary(fit.npb.14)
npb.sum.14$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.2208720 6.711495 -22.803458 1.487951 0.234
## [2,] -21.1917922 20.300690 -64.503294 0.000000 0.748
## [3,] -1.5420847 6.081411 -22.659259 2.424170 0.200
## [4,] -3.3387828 9.164850 -32.417900 0.000000 0.254
## [5,] -1.2164232 6.085029 -18.771920 6.704256 0.202
## [6,] -0.3928061 4.866547 -12.933445 8.727297 0.172
## [7,] -0.1886759 4.562750 -10.614897 9.006422 0.186
## [8,] 1.9671743 9.803766 -9.345317 38.040782 0.210
## [9,] 0.3059293 5.512448 -10.409413 13.616308 0.188
## [10,] -1.0485428 5.747603 -16.448024 6.972009 0.210
## [11,] -4.8200773 10.904677 -38.663265 3.040430 0.326
## [12,] -2.2364080 8.259488 -24.910237 5.348292 0.268
## [13,] -5.2746140 12.231653 -42.316979 1.903118 0.320
## [14,] -0.9090102 4.871460 -13.365501 2.922642 0.174
## [15,] -3.5391613 8.769814 -31.460777 2.242225 0.292
## [16,] -3.4732758 9.609786 -33.340210 5.511336 0.310
## [17,] -10.7033356 17.707997 -60.140820 0.000000 0.446
## [18,] -2.3657150 7.688985 -26.275109 3.211347 0.240
## [19,] -2.4531044 8.486316 -25.800574 5.974371 0.238plot(fit.npb.14$beta[,1], type = "l")
plot(fit.npb.14$beta[,2], type = "l")
plot(fit.npb.14$beta[,13], type = "l")
4.1.4.4 Try making alpha.pi and beta.pi 10
This might be getting too big… Is there enough space in the tails?
plot(density(rbeta(10000, 10, 10)))
priors.npb.15 <- list(alpha.pi = 10, beta.pi = 10, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 1)
fit.npb.15 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.15, interact = F)
npb.sum.15 <- summary(fit.npb.15)
npb.sum.15$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.8857311 7.232977 -23.56255 5.3865842 0.306
## [2,] -17.1289824 18.507669 -63.63778 0.0000000 0.712
## [3,] -1.6323034 6.176281 -16.91383 6.2286586 0.238
## [4,] -3.6020419 8.267143 -29.03935 2.8766411 0.332
## [5,] -1.5654875 5.200797 -16.56090 3.6318501 0.214
## [6,] -0.6052462 7.281443 -17.65679 16.2820039 0.224
## [7,] -0.2546045 5.397410 -13.79061 12.1249885 0.198
## [8,] 1.5759810 9.269974 -12.00373 37.5717956 0.224
## [9,] -0.3945213 6.360807 -14.74431 11.4964058 0.214
## [10,] -1.2363866 5.434546 -16.68779 5.3865842 0.212
## [11,] -5.1448347 9.714039 -31.62555 0.8892998 0.362
## [12,] -3.1353175 7.871987 -26.39008 4.2840909 0.296
## [13,] -4.7910311 9.203093 -29.48755 0.7573745 0.382
## [14,] -0.9290989 6.399126 -14.91728 9.2503359 0.246
## [15,] -3.1021597 7.106243 -23.80494 0.6839960 0.288
## [16,] -4.2095545 8.910515 -30.64462 3.1182682 0.354
## [17,] -9.8871161 14.791150 -51.45763 0.0000000 0.516
## [18,] -2.5865340 7.240974 -25.16143 2.8766411 0.288
## [19,] -2.6717984 8.003829 -24.29969 3.6711970 0.280plot(fit.npb.15$beta[,1], type = "l")
plot(fit.npb.15$beta[,2], type = "l")
plot(fit.npb.15$beta[,13], type = "l")
4.1.5 Set alpha.pi and beta.pi to 8, readjust a.phi1 and sig2inv.mu1
Set alpha.pi and beta.pi to 8, try adjusting a.phi1 and sig2inv.mu1
4.1.5.1 Try making a.phi1 = 10 and sig2inv.mu1 = 1
priors.npb.16 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 1)
fit.npb.16 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.16, interact = F)
npb.sum.16 <- summary(fit.npb.16)
npb.sum.16$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.5584270 7.366529 -23.113622 3.844689 0.274
## [2,] -25.6541530 20.781097 -70.525111 0.000000 0.822
## [3,] -0.8456003 5.351849 -16.155950 6.628430 0.186
## [4,] -3.9716715 9.687439 -33.874152 3.100136 0.316
## [5,] -0.6566620 5.421212 -15.942915 6.547565 0.202
## [6,] -0.5472398 5.943481 -13.675801 12.959635 0.222
## [7,] 0.5289133 7.935148 -12.285135 19.441791 0.220
## [8,] 3.6669225 11.767144 -6.936092 44.027025 0.250
## [9,] 0.9840048 7.205198 -11.178864 22.955402 0.224
## [10,] -0.8924849 6.122219 -16.471644 9.911731 0.208
## [11,] -5.2796511 11.725045 -39.869035 2.702883 0.346
## [12,] -2.9070741 9.647720 -31.198810 8.953711 0.298
## [13,] -5.6108588 12.693757 -42.736833 6.458987 0.360
## [14,] 0.1968864 9.101304 -16.269186 21.797660 0.232
## [15,] -3.9612949 10.534515 -38.442630 3.788611 0.298
## [16,] -2.7523511 8.514096 -28.261148 5.919833 0.280
## [17,] -13.8663411 18.848322 -57.570802 0.123517 0.534
## [18,] -1.3469540 6.406979 -17.682461 9.557808 0.250
## [19,] -2.8201205 9.091258 -30.279149 5.574324 0.270plot(fit.npb.16$beta[,1], type = "l")
plot(fit.npb.16$beta[,2], type = "l")
plot(fit.npb.16$beta[,13], type = "l")
Same priors, but now with scaleY = F
priors.npb.16a <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 1)
fit.npb.16a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = FALSE,
priors = priors.npb.16a, interact = F)
npb.sum.16a <- summary(fit.npb.16a)
npb.sum.16a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.01991776141 0.2863296 -0.6992964 0.5727729 0.464
## [2,] 0.00056061079 0.3187154 -0.7958305 0.7398862 0.506
## [3,] -0.03471161517 0.3055429 -0.7110587 0.6485434 0.480
## [4,] -0.02622573229 0.2940953 -0.6886811 0.6722611 0.494
## [5,] -0.03089337315 0.3053399 -0.7938389 0.6140577 0.486
## [6,] -0.01576318329 0.3080216 -0.7080183 0.6828620 0.494
## [7,] -0.03381090913 0.3125751 -0.8258076 0.6287432 0.478
## [8,] -0.00861707349 0.3123751 -0.7661726 0.6777252 0.518
## [9,] -0.02157509124 0.3382525 -0.8302362 0.6707648 0.546
## [10,] -0.02870830928 0.3003430 -0.7902254 0.6221704 0.524
## [11,] -0.02532884091 0.3069601 -0.7470368 0.6530465 0.502
## [12,] -0.02095698714 0.3255910 -0.7847839 0.6771903 0.498
## [13,] -0.02059436427 0.2980543 -0.7584723 0.6586088 0.502
## [14,] -0.01576474199 0.3142070 -0.6979910 0.6499558 0.510
## [15,] -0.02125723989 0.3174181 -0.7684473 0.6804516 0.488
## [16,] 0.00006688786 0.3072310 -0.7088393 0.7004202 0.486
## [17,] -0.02806268605 0.3062034 -0.7106627 0.6624262 0.488
## [18,] -0.00615786566 0.3145777 -0.7314019 0.6625514 0.492
## [19,] -0.02278172580 0.3291935 -0.7219907 0.7325880 0.512plot(fit.npb.16a$beta[,1], type = "l")
plot(fit.npb.16a$beta[,2], type = "l")
plot(fit.npb.16a$beta[,13], type = "l")
4.1.5.2 Try making a.phi1 = 10 and sig2inv.mu1 = 10
priors.npb.17 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 10)
fit.npb.17 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.17, interact = F)
npb.sum.17 <- summary(fit.npb.17)
npb.sum.17$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.1114031 7.367810 -25.931299 9.363107 0.274
## [2,] -25.3346069 19.863539 -66.334946 0.000000 0.848
## [3,] -1.2404926 5.768699 -17.263289 9.161327 0.256
## [4,] -3.8940779 9.888461 -32.943501 5.677994 0.312
## [5,] -0.6420384 5.274159 -13.895186 9.204381 0.230
## [6,] -0.6245669 6.644525 -17.225677 13.345719 0.256
## [7,] 0.3270219 7.875332 -13.175990 17.337589 0.250
## [8,] 4.2127861 13.452337 -10.933251 47.477502 0.306
## [9,] 1.1309278 7.689213 -9.963829 23.003936 0.242
## [10,] -1.5547519 7.118842 -19.923356 11.966823 0.242
## [11,] -6.8228893 13.060407 -42.538875 5.610858 0.432
## [12,] -2.2978183 9.767088 -28.140006 11.749363 0.304
## [13,] -5.9216198 12.252662 -41.722235 6.660095 0.394
## [14,] -0.3628362 6.824363 -15.578790 15.121946 0.240
## [15,] -3.8462739 9.954491 -31.914553 6.102279 0.336
## [16,] -3.0798149 9.305046 -29.956764 10.470342 0.340
## [17,] -14.8435061 17.959995 -56.317799 0.000000 0.592
## [18,] -1.5879665 7.736488 -21.215828 11.581411 0.268
## [19,] -2.7080880 8.515935 -27.715093 8.802932 0.302plot(fit.npb.17$beta[,1], type = "l")
plot(fit.npb.17$beta[,2], type = "l")
plot(fit.npb.17$beta[,13], type = "l")
4.1.5.3 Try making a.phi1 = 10 and sig2inv.mu1 = 100
priors.npb.18 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 100)
fit.npb.18 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.18, interact = F)
npb.sum.18 <- summary(fit.npb.18)
npb.sum.18$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.7940023 7.665154 -25.139592 1.967637 0.262
## [2,] -24.6794003 19.177369 -64.105042 0.000000 0.826
## [3,] -2.1369699 7.197838 -23.287142 2.457191 0.234
## [4,] -4.0511400 9.299041 -31.811440 2.462465 0.330
## [5,] -0.6707083 5.768537 -16.919951 7.073050 0.168
## [6,] -0.6148905 6.696902 -17.526985 15.081068 0.212
## [7,] -0.1515304 8.041328 -15.751521 17.158266 0.226
## [8,] 5.4699606 14.981968 -5.774078 51.032414 0.246
## [9,] 0.6457981 7.848650 -11.650429 30.095922 0.180
## [10,] -1.4251170 7.085737 -22.208937 7.116970 0.228
## [11,] -6.3271767 12.483270 -43.419739 2.008763 0.380
## [12,] -3.3099004 10.237558 -32.374922 4.621957 0.266
## [13,] -5.9102346 11.934866 -39.143598 1.522695 0.350
## [14,] 0.1825934 11.150450 -19.088496 31.667661 0.236
## [15,] -4.5956090 10.565888 -30.641680 1.299086 0.316
## [16,] -3.7256016 10.659967 -32.794314 5.974629 0.308
## [17,] -12.7496421 17.777733 -59.090938 0.000000 0.526
## [18,] -2.1534587 7.856944 -25.086513 3.481331 0.248
## [19,] -2.0881613 8.620420 -24.424391 3.983888 0.222plot(fit.npb.18$beta[,1], type = "l")
plot(fit.npb.18$beta[,2], type = "l")
plot(fit.npb.18$beta[,13], type = "l")
4.1.5.4 Try making a.phi1 = 100 and sig2inv.mu1 = 1
priors.npb.19 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 1)
fit.npb.19 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.19, interact = F)
npb.sum.19 <- summary(fit.npb.19)
npb.sum.19$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.93945571 7.436919 -24.014768 11.1380086 0.290
## [2,] -29.35356001 20.315233 -69.971480 0.0000000 0.886
## [3,] -0.92912783 8.390025 -23.927612 17.8853035 0.290
## [4,] -4.24255171 11.019911 -37.013982 11.3654517 0.352
## [5,] -0.02346836 7.140468 -17.342551 19.0096949 0.248
## [6,] 0.25407353 7.350061 -14.627499 19.5678783 0.240
## [7,] 1.57211236 8.970933 -14.205714 27.7716715 0.280
## [8,] 5.84507537 14.003653 -9.061805 47.5658242 0.364
## [9,] 3.04650766 10.489651 -10.050578 37.3044036 0.302
## [10,] -0.62670622 8.574075 -20.521278 20.0444767 0.290
## [11,] -6.73510501 13.217648 -44.885409 7.8525053 0.430
## [12,] -2.04120535 9.872410 -29.555056 17.1520514 0.326
## [13,] -6.11778167 12.998699 -40.797392 7.4243984 0.386
## [14,] 0.93850331 8.248327 -12.573063 24.4633603 0.286
## [15,] -4.23998411 10.449073 -33.210287 9.0480373 0.374
## [16,] -3.66326995 10.544825 -34.165383 10.0206015 0.338
## [17,] -14.96004205 18.360076 -58.461572 0.4857887 0.612
## [18,] -1.47046601 10.311921 -26.579843 20.5337795 0.336
## [19,] -2.68239003 10.918269 -33.565613 17.0679189 0.334plot(fit.npb.19$beta[,1], type = "l")
plot(fit.npb.19$beta[,2], type = "l")
plot(fit.npb.19$beta[,13], type = "l")
4.1.5.5 Try making a.phi1 = 100 and sig2inv.mu1 = 10
priors.npb.20 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 10)
fit.npb.20 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.20, interact = F)
npb.sum.20 <- summary(fit.npb.20)
npb.sum.20$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.0957585 7.253974 -19.996235 12.133194 0.260
## [2,] -30.7165080 20.599113 -70.129391 0.000000 0.876
## [3,] -1.1419660 8.278356 -23.608839 17.259197 0.272
## [4,] -2.8683240 7.995775 -28.839509 4.979373 0.284
## [5,] -0.2349687 7.355764 -15.296693 19.075479 0.248
## [6,] 0.4651001 7.262474 -14.319838 21.345014 0.268
## [7,] 2.3831534 9.479152 -11.090020 29.683408 0.270
## [8,] 6.7562698 14.483801 -8.564801 48.282371 0.388
## [9,] 1.3374750 8.602574 -11.469325 30.417517 0.262
## [10,] -0.4535202 7.685077 -20.453290 19.872858 0.264
## [11,] -6.8788605 13.754454 -45.639225 5.972641 0.412
## [12,] -1.2988578 8.958916 -23.994833 13.315457 0.296
## [13,] -6.4462750 13.404348 -43.699827 8.962873 0.420
## [14,] 0.4766060 7.309769 -14.083425 20.419289 0.254
## [15,] -4.1468831 10.703679 -34.713864 10.198616 0.374
## [16,] -3.7334116 11.038755 -35.748628 12.120982 0.338
## [17,] -14.5536692 19.178035 -59.659277 0.000000 0.556
## [18,] -1.4283640 8.830606 -24.765649 16.967138 0.288
## [19,] -2.6752086 11.080455 -39.129165 16.421379 0.330plot(fit.npb.20$beta[,1], type = "l")
plot(fit.npb.20$beta[,2], type = "l")
plot(fit.npb.20$beta[,13], type = "l")
4.1.5.6 Try making a.phi1 = 100 and sig2inv.mu1 = 100
priors.npb.21 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 100)
fit.npb.21 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.21, interact = F)
npb.sum.21 <- summary(fit.npb.21)
npb.sum.21$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.598852431 8.120628 -22.125745 14.6291093 0.350
## [2,] -28.636717497 20.706557 -70.367184 0.0000000 0.894
## [3,] -0.615990793 8.153533 -18.892649 16.9367377 0.322
## [4,] -4.058180165 9.708837 -33.927978 10.4658336 0.426
## [5,] -0.003474108 7.731554 -15.980152 23.1872641 0.322
## [6,] 0.551072649 8.516808 -16.233849 23.5174997 0.296
## [7,] 1.834206486 9.568153 -14.109403 29.9184057 0.338
## [8,] 6.766290708 14.955496 -8.186632 47.3844319 0.414
## [9,] 1.751247526 8.799965 -10.038815 26.8865154 0.290
## [10,] -1.519519066 7.488465 -21.309451 14.6719829 0.318
## [11,] -5.018906137 11.662643 -37.760597 11.9308543 0.434
## [12,] -1.677747508 11.211820 -27.835487 21.5543879 0.394
## [13,] -4.838752463 11.193045 -36.043497 10.4594162 0.430
## [14,] 0.863907073 8.691483 -13.560335 26.9581192 0.328
## [15,] -4.327300531 9.899823 -31.088168 5.3204325 0.386
## [16,] -3.276965016 10.419959 -31.561708 16.7571722 0.424
## [17,] -14.099124604 18.520190 -62.348545 0.1831654 0.598
## [18,] -0.759890932 8.582973 -18.788232 20.9963292 0.334
## [19,] -3.076637438 9.461079 -31.613003 12.5061512 0.360plot(fit.npb.21$beta[,1], type = "l")
plot(fit.npb.21$beta[,2], type = "l")
plot(fit.npb.21$beta[,13], type = "l")
Same priors, but with scaleY = F
priors.npb.21a <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 100)
fit.npb.21a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = F,
priors = priors.npb.21a, interact = F)
npb.sum.21a <- summary(fit.npb.21a)
npb.sum.21a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] 0.0051412432 0.08469866 -0.1828181 0.2085003 0.502
## [2,] 0.0037644320 0.09304471 -0.2130544 0.2124937 0.540
## [3,] 0.0026560704 0.08663993 -0.1781822 0.2117353 0.530
## [4,] 0.0013612505 0.07993577 -0.1711556 0.1806295 0.502
## [5,] -0.0023618224 0.08429401 -0.2003212 0.1807728 0.550
## [6,] 0.0028257961 0.07825187 -0.1775277 0.2066858 0.486
## [7,] -0.0025310235 0.08453026 -0.1876556 0.1875710 0.538
## [8,] 0.0013732366 0.08632333 -0.1878732 0.1925771 0.518
## [9,] 0.0033625563 0.08540797 -0.1894751 0.1967161 0.520
## [10,] -0.0049627231 0.08504732 -0.2075164 0.1905664 0.556
## [11,] -0.0038248546 0.08145215 -0.1857765 0.1905289 0.514
## [12,] -0.0073032570 0.08269624 -0.1918137 0.1881970 0.510
## [13,] -0.0037368492 0.08167179 -0.1870376 0.1790697 0.504
## [14,] 0.0046757747 0.08909916 -0.1864469 0.2050904 0.574
## [15,] -0.0021973315 0.08784817 -0.1918137 0.1872830 0.512
## [16,] 0.0019503104 0.08383080 -0.1736276 0.2080948 0.530
## [17,] 0.0077824576 0.08742294 -0.1830547 0.2071179 0.540
## [18,] 0.0001896831 0.08418336 -0.1850637 0.1905664 0.524
## [19,] 0.0063805729 0.08233081 -0.1782442 0.2015107 0.492plot(fit.npb.21a$beta[,1], type = "l")
plot(fit.npb.21a$beta[,2], type = "l")
plot(fit.npb.21a$beta[,13], type = "l")
4.1.5.7 Try making a.phi1 = 10e5 and sig2inv.mu1 = 10e5
priors.npb.22 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10e5, sig2inv.mu1 = 10e5)
fit.npb.22 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.22, interact = F)
npb.sum.22 <- summary(fit.npb.22)
npb.sum.22$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.035306284 0.4087304 -1.0444576 0.8734612 0.520
## [2,] -0.047707128 0.4067405 -0.9722595 0.8787816 0.502
## [3,] -0.035577082 0.4012861 -1.0507546 0.8651339 0.506
## [4,] -0.057457586 0.4065833 -1.0596305 0.7537023 0.526
## [5,] -0.044657655 0.4247126 -1.1437954 0.9025106 0.494
## [6,] -0.047098218 0.4545910 -1.1189796 1.0003863 0.506
## [7,] -0.002177765 0.4795466 -1.0556023 1.0804129 0.558
## [8,] -0.011803902 0.3942380 -0.9792511 0.9517408 0.498
## [9,] -0.033054139 0.4535246 -1.1286721 0.9043937 0.548
## [10,] -0.043970919 0.4285394 -1.1832949 0.9432537 0.510
## [11,] -0.034836820 0.4687112 -1.0535954 1.1372874 0.526
## [12,] -0.055681433 0.4141023 -1.1090113 0.8180197 0.478
## [13,] -0.061520085 0.4261508 -1.1048181 0.8484011 0.506
## [14,] -0.014160757 0.4547002 -1.0934082 0.9126038 0.534
## [15,] -0.047145647 0.4326528 -1.1168833 0.8238743 0.478
## [16,] -0.048599028 0.4232817 -1.0691288 0.7790946 0.528
## [17,] -0.070057692 0.4356459 -1.1557702 0.7609288 0.536
## [18,] -0.020791184 0.4355400 -1.1591367 0.9563392 0.476
## [19,] -0.040411595 0.4550934 -1.0530173 0.9865404 0.506plot(fit.npb.22$beta[,1], type = "l")
plot(fit.npb.22$beta[,2], type = "l")
plot(fit.npb.22$beta[,13], type = "l")
Same priors, but with scaleY = F
priors.npb.22a <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10e5, sig2inv.mu1 = 10e5)
fit.npb.22a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = FALSE,
priors = priors.npb.22a, interact = F)
npb.sum.22a <- summary(fit.npb.22a)
npb.sum.22a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] 0.000006182868 0.0008619421 -0.001979458 0.002071416 0.520
## [2,] -0.000052064014 0.0008040348 -0.001874284 0.001670370 0.524
## [3,] -0.000046371401 0.0008962762 -0.002141472 0.001883780 0.522
## [4,] 0.000039298564 0.0009298157 -0.002246873 0.001975836 0.522
## [5,] -0.000017857840 0.0008660585 -0.002082111 0.001861066 0.546
## [6,] -0.000076962013 0.0009002153 -0.002196562 0.001778247 0.512
## [7,] -0.000017274586 0.0008623192 -0.001970859 0.002013509 0.528
## [8,] -0.000042419780 0.0008878686 -0.001940788 0.001864174 0.528
## [9,] -0.000037825550 0.0008999260 -0.002129650 0.001889525 0.542
## [10,] -0.000071637744 0.0008846355 -0.002220330 0.001768216 0.506
## [11,] 0.000018945719 0.0009117751 -0.002224579 0.002037273 0.520
## [12,] 0.000076549490 0.0007989603 -0.002009522 0.001835409 0.470
## [13,] 0.000065717471 0.0008749839 -0.001984107 0.001916398 0.524
## [14,] -0.000030690991 0.0008669480 -0.002107817 0.001895524 0.550
## [15,] 0.000019370894 0.0009155382 -0.001984431 0.001988277 0.542
## [16,] -0.000036712651 0.0009179785 -0.002146234 0.001986838 0.502
## [17,] -0.000028941612 0.0008525617 -0.001949226 0.001762325 0.512
## [18,] 0.000010433393 0.0009188860 -0.002237139 0.002022782 0.514
## [19,] -0.000023885689 0.0008964864 -0.002187121 0.001901227 0.540plot(fit.npb.22a$beta[,1], type = "l")
plot(fit.npb.22a$beta[,2], type = "l")
plot(fit.npb.22a$beta[,13], type = "l")
4.1.6 Set alpha.pi and beta.pi to 5, readjust a.phi1 and sig2inv.mu1
Set alpha.pi and beta.pi to 5, rather than 8, and try adjusting a.phi1 and sig2inv.mu1
4.1.6.1 Try making a.phi1 = 10 and sig2inv.mu1 = 1
priors.npb.23 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 1)
fit.npb.23 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.23, interact = F)
npb.sum.23 <- summary(fit.npb.23)
npb.sum.23$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.2669275 6.818737 -23.985602 7.226331 0.284
## [2,] -23.8968298 21.256605 -70.757009 0.000000 0.802
## [3,] -1.3764332 5.783617 -18.251351 6.223349 0.220
## [4,] -3.1824556 8.462735 -27.548227 3.788901 0.308
## [5,] -1.1498466 5.614243 -17.671669 6.601940 0.232
## [6,] -1.0777643 6.557314 -19.429210 11.728310 0.250
## [7,] 0.3905512 8.642639 -14.259111 24.121953 0.232
## [8,] 3.6174570 12.688717 -9.551844 45.028736 0.256
## [9,] 0.8883675 8.857065 -13.113582 23.017336 0.236
## [10,] -0.5558551 5.804998 -14.312065 8.953615 0.212
## [11,] -4.8874443 10.554609 -35.472399 5.931729 0.372
## [12,] -2.0602220 7.428268 -24.099109 8.386962 0.280
## [13,] -5.3142433 10.814921 -40.185742 1.913820 0.358
## [14,] -0.6998181 6.073537 -16.303471 11.832843 0.236
## [15,] -2.9213200 7.868612 -26.124243 4.919153 0.308
## [16,] -3.8903774 10.524390 -35.856330 8.674976 0.354
## [17,] -11.2928049 17.090541 -56.513393 1.655472 0.514
## [18,] -1.8947305 7.612136 -22.642968 8.646449 0.276
## [19,] -3.2289844 8.910317 -28.874166 3.991413 0.294plot(fit.npb.23$beta[,1], type = "l")
plot(fit.npb.23$beta[,2], type = "l")
plot(fit.npb.23$beta[,13], type = "l")
Same priors, but now with scaleY = F
priors.npb.23a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 1)
fit.npb.23a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = FALSE,
priors = priors.npb.23a, interact = F)
npb.sum.23a <- summary(fit.npb.23a)
npb.sum.23a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.0007196103 0.2964622 -0.6553441 0.7020867 0.510
## [2,] -0.0120527174 0.2632085 -0.6428513 0.5712860 0.498
## [3,] -0.0003302213 0.2964493 -0.6474469 0.6524802 0.538
## [4,] 0.0038542799 0.2704418 -0.6072693 0.5571160 0.502
## [5,] 0.0073328019 0.2603254 -0.5727369 0.5508989 0.472
## [6,] 0.0149783637 0.2940539 -0.6596654 0.7206298 0.482
## [7,] -0.0148465978 0.3064579 -0.6605561 0.6640202 0.534
## [8,] 0.0042084310 0.3055809 -0.6986858 0.6806700 0.502
## [9,] 0.0040830043 0.2701805 -0.5430343 0.6331375 0.480
## [10,] 0.0070490226 0.2993666 -0.6461819 0.7809095 0.516
## [11,] 0.0104681999 0.2989888 -0.6058662 0.7507703 0.524
## [12,] -0.0009639157 0.2873621 -0.6392656 0.6505214 0.512
## [13,] -0.0202041392 0.2846702 -0.6886220 0.5432586 0.532
## [14,] 0.0097133647 0.2979735 -0.6236524 0.7223282 0.506
## [15,] 0.0268837151 0.3114360 -0.6323325 0.7621870 0.508
## [16,] -0.0019335793 0.2594860 -0.6081490 0.6522747 0.462
## [17,] -0.0144824164 0.2984423 -0.7850978 0.6364156 0.514
## [18,] 0.0088153563 0.2645159 -0.5405364 0.7058493 0.468
## [19,] 0.0197624609 0.3104259 -0.6356867 0.7961824 0.490plot(fit.npb.23a$beta[,1], type = "l")
plot(fit.npb.23a$beta[,2], type = "l")
plot(fit.npb.23a$beta[,13], type = "l")
4.1.6.2 Try making a.phi1 = 10 and sig2inv.mu1 = 10
priors.npb.24 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 10)
fit.npb.24 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.24, interact = F)
npb.sum.24 <- summary(fit.npb.24)
npb.sum.24$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.1744090 6.356825 -21.365240 6.5471919 0.180
## [2,] -25.9222430 21.582287 -70.704025 0.0000000 0.764
## [3,] -0.8719158 5.355631 -17.380790 5.0388036 0.154
## [4,] -2.6852993 8.750953 -31.977107 0.8560601 0.214
## [5,] -0.3376489 4.707140 -13.656701 4.3798985 0.122
## [6,] -0.1516117 6.011091 -12.307358 11.9384392 0.158
## [7,] 0.2033702 5.280983 -9.111321 13.6905185 0.150
## [8,] 4.1573659 12.676688 -5.377147 45.2223780 0.212
## [9,] 0.8289965 7.574134 -10.641479 24.7475649 0.164
## [10,] -0.6102648 6.032846 -16.157754 10.9900904 0.156
## [11,] -4.5080221 10.836041 -37.483169 2.7381582 0.274
## [12,] -1.5336841 7.601180 -23.358433 3.3745605 0.192
## [13,] -4.5377605 11.415565 -39.490858 0.9601621 0.280
## [14,] -0.1297684 6.641396 -12.858573 18.9852703 0.178
## [15,] -3.8896588 10.341413 -37.505833 2.9042601 0.282
## [16,] -3.6972612 9.804580 -36.883240 1.6065287 0.230
## [17,] -11.4272777 18.220939 -58.972882 0.0000000 0.434
## [18,] -1.5892077 7.539632 -22.211002 8.4387392 0.204
## [19,] -2.7649340 9.661202 -36.632517 3.6407285 0.210plot(fit.npb.24$beta[,1], type = "l")
plot(fit.npb.24$beta[,2], type = "l")
plot(fit.npb.24$beta[,13], type = "l")
4.1.6.3 Try making a.phi1 = 10 and sig2inv.mu1 = 100
priors.npb.25 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10, sig2inv.mu1 = 100)
fit.npb.25 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.25, interact = F)
npb.sum.25 <- summary(fit.npb.25)
npb.sum.25$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.66648315 7.226089 -26.871111 6.1849869 0.200
## [2,] -28.80748355 20.864602 -67.769965 0.0000000 0.838
## [3,] -1.42558012 6.522784 -22.771262 4.4281229 0.192
## [4,] -2.76921892 8.454176 -31.199743 1.5756317 0.216
## [5,] -0.28525201 4.552973 -10.345389 7.3140781 0.144
## [6,] -0.02273501 5.207522 -10.280686 12.2972152 0.128
## [7,] 0.42838036 6.509331 -8.832926 15.5831741 0.146
## [8,] 3.61899464 11.611365 -4.922887 41.1492615 0.204
## [9,] 0.65651537 6.752541 -9.044913 17.4790875 0.154
## [10,] -0.28000249 4.486445 -11.157330 7.2672729 0.164
## [11,] -5.83525887 12.531493 -42.110405 2.3667362 0.306
## [12,] -2.20787228 8.966495 -29.191497 3.2932711 0.220
## [13,] -7.12779870 14.099981 -48.626023 1.9446005 0.344
## [14,] 0.36657831 8.642415 -15.834819 24.3479811 0.194
## [15,] -4.05853291 10.317329 -35.917289 0.7062024 0.242
## [16,] -3.19567962 9.827098 -34.041145 4.7932942 0.244
## [17,] -12.25439359 18.517567 -58.794296 0.0000000 0.434
## [18,] -1.90909272 8.525734 -27.250940 5.5485470 0.204
## [19,] -3.46970786 10.432867 -37.223133 3.3684512 0.248plot(fit.npb.25$beta[,1], type = "l")
plot(fit.npb.25$beta[,2], type = "l")
plot(fit.npb.25$beta[,13], type = "l")
4.1.6.4 Try making a.phi1 = 100 and sig2inv.mu1 = 1
priors.npb.26 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 1)
fit.npb.26 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.26, interact = F)
npb.sum.26 <- summary(fit.npb.26)
npb.sum.26$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.2621545 7.405051 -22.005981 11.964441 0.262
## [2,] -31.8730363 20.471567 -69.939915 0.000000 0.862
## [3,] -0.7504150 8.021987 -20.950051 16.075527 0.282
## [4,] -4.1027594 11.475524 -41.354956 9.531501 0.336
## [5,] 0.1652421 5.658890 -14.033363 14.483537 0.236
## [6,] 0.7633474 7.314781 -14.251274 18.801390 0.296
## [7,] 2.1124896 9.166049 -11.404178 28.577810 0.266
## [8,] 7.6878235 15.418804 -7.034124 50.178836 0.410
## [9,] 2.7160023 9.653674 -9.890912 31.315221 0.314
## [10,] 0.3323067 8.778377 -17.837627 18.701679 0.288
## [11,] -6.5489320 13.681433 -44.351904 6.346462 0.378
## [12,] -2.2381381 9.969420 -30.264782 15.145125 0.320
## [13,] -4.4238949 11.719498 -37.106223 9.696419 0.354
## [14,] 1.1390514 7.813899 -13.346768 22.660832 0.264
## [15,] -3.6777909 10.307871 -34.607732 9.943842 0.330
## [16,] -1.4154880 10.777865 -31.506473 19.824269 0.322
## [17,] -16.4797330 20.840367 -63.168555 1.184975 0.564
## [18,] -0.6581118 8.235614 -23.385310 18.644877 0.318
## [19,] -1.6463141 9.357480 -30.512594 13.744361 0.288plot(fit.npb.26$beta[,1], type = "l")
plot(fit.npb.26$beta[,2], type = "l")
plot(fit.npb.26$beta[,13], type = "l")
4.1.6.5 Try making a.phi1 = 100 and sig2inv.mu1 = 10
priors.npb.27 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 10)
fit.npb.27 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.27, interact = F)
npb.sum.27 <- summary(fit.npb.27)
npb.sum.27$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.01163147 6.974847 -23.305300 12.063268 0.248
## [2,] -30.66127687 20.760120 -70.934734 0.000000 0.872
## [3,] -1.00436024 7.982394 -24.827645 15.695069 0.270
## [4,] -3.22484204 10.981124 -35.503109 10.899975 0.304
## [5,] 0.65483898 7.215895 -14.305391 21.422911 0.226
## [6,] 1.05202716 8.610743 -13.386552 26.288227 0.240
## [7,] 1.04745993 9.161936 -18.039125 24.969603 0.270
## [8,] 7.61348945 15.349643 -9.331034 48.223341 0.388
## [9,] 2.22017934 9.675796 -14.805064 31.309280 0.256
## [10,] 0.06297748 7.844041 -17.031138 18.704314 0.260
## [11,] -6.49365289 13.166477 -44.745196 5.150326 0.368
## [12,] -1.23067543 10.029990 -29.537828 21.137705 0.272
## [13,] -5.58891305 13.398458 -44.759448 6.707673 0.390
## [14,] 0.48784423 9.658446 -21.374144 27.318140 0.246
## [15,] -3.70084647 9.962155 -34.880874 3.908209 0.288
## [16,] -1.89826339 10.065355 -26.207320 15.512113 0.304
## [17,] -14.86352388 19.333278 -60.535336 0.000000 0.540
## [18,] -1.82229314 8.936614 -29.039057 12.853870 0.270
## [19,] -3.11206478 11.992936 -39.216810 15.476387 0.326plot(fit.npb.27$beta[,1], type = "l")
plot(fit.npb.27$beta[,2], type = "l")
plot(fit.npb.27$beta[,13], type = "l")
4.1.6.6 Try making a.phi1 = 100 and sig2inv.mu1 = 100
priors.npb.28 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 100)
fit.npb.28 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.28, interact = F)
npb.sum.28 <- summary(fit.npb.28)
npb.sum.28$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.5967933 7.684725 -23.845516 12.9857927 0.262
## [2,] -30.3509594 21.129884 -69.873352 0.0000000 0.848
## [3,] -1.1081974 7.846097 -22.418736 14.0079554 0.260
## [4,] -3.1540371 9.566356 -32.632348 8.9516887 0.292
## [5,] -0.1257847 6.959253 -16.825701 18.5528850 0.254
## [6,] 0.5449194 7.866129 -13.596188 22.6291587 0.288
## [7,] 1.1863816 9.553280 -13.924976 27.4139131 0.288
## [8,] 8.7649384 17.034173 -10.170230 53.1274332 0.444
## [9,] 3.5466042 11.407061 -9.182936 41.9792743 0.328
## [10,] -0.3097942 7.638762 -19.690696 17.8438100 0.274
## [11,] -5.9581373 12.160559 -40.851977 6.9929063 0.382
## [12,] -1.6259110 10.051129 -26.222054 19.4293422 0.316
## [13,] -6.1755473 13.010573 -42.585762 7.8574922 0.394
## [14,] 1.7462565 10.147823 -17.436863 33.7645626 0.272
## [15,] -3.6360794 10.305510 -36.022150 9.1708419 0.324
## [16,] -1.4794849 11.214146 -30.263661 21.3994828 0.330
## [17,] -17.2136283 20.264874 -63.256164 0.4187659 0.618
## [18,] -1.0341528 9.680756 -27.793811 19.1325054 0.306
## [19,] -2.1742274 11.416120 -33.177502 19.8048518 0.328plot(fit.npb.28$beta[,1], type = "l")
plot(fit.npb.28$beta[,2], type = "l")
plot(fit.npb.28$beta[,13], type = "l")
Same priors, but with scaleY = F
priors.npb.28a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 100, sig2inv.mu1 = 100)
fit.npb.28a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = F,
priors = priors.npb.28a, interact = F)
npb.sum.28a <- summary(fit.npb.28a)
npb.sum.28a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] 0.00304819944 0.07585417 -0.1791721 0.1588592 0.504
## [2,] 0.00368920822 0.08356244 -0.1832763 0.1953595 0.484
## [3,] 0.00230653550 0.08395801 -0.1981486 0.1913326 0.518
## [4,] 0.00349433034 0.08395506 -0.1834283 0.2139654 0.484
## [5,] -0.00392341030 0.07385651 -0.1812940 0.1651630 0.480
## [6,] -0.00008981685 0.08121901 -0.1920336 0.1843800 0.484
## [7,] 0.00700464294 0.08178275 -0.1759053 0.2133049 0.480
## [8,] 0.00157109764 0.07703858 -0.1735443 0.1754262 0.476
## [9,] 0.01105053383 0.07811738 -0.1558621 0.2026439 0.496
## [10,] 0.00674268321 0.08007177 -0.1618587 0.1982301 0.484
## [11,] 0.00070700654 0.07406451 -0.1659681 0.1679017 0.482
## [12,] 0.00148550141 0.08263396 -0.1969318 0.1810100 0.502
## [13,] 0.00514617985 0.08217526 -0.1626479 0.1919755 0.514
## [14,] -0.00020879633 0.08166007 -0.1723812 0.1923978 0.502
## [15,] -0.00140959112 0.07936366 -0.1735443 0.1720113 0.488
## [16,] 0.00212519609 0.07457426 -0.1619751 0.1777630 0.494
## [17,] 0.00486790268 0.07542233 -0.1648815 0.1680491 0.516
## [18,] -0.00225180822 0.08085503 -0.1661710 0.1791710 0.516
## [19,] 0.00332850383 0.07821127 -0.1611759 0.1848179 0.520plot(fit.npb.28a$beta[,1], type = "l")
plot(fit.npb.28a$beta[,2], type = "l")
plot(fit.npb.28a$beta[,13], type = "l")
4.1.6.7 Try making a.phi1 = 10e5 and sig2inv.mu1 = 10e5
priors.npb.29 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10e5, sig2inv.mu1 = 10e5)
fit.npb.29 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.29, interact = F)
npb.sum.29 <- summary(fit.npb.29)
npb.sum.29$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.006677226 0.4192719 -0.9510673 0.9945879 0.492
## [2,] -0.069118614 0.4012581 -1.1079014 0.7428316 0.484
## [3,] -0.038108583 0.4021976 -1.0728872 0.8241823 0.506
## [4,] -0.048976955 0.4261434 -1.0467512 0.8276469 0.532
## [5,] -0.053049503 0.4064875 -1.0925019 0.7353015 0.498
## [6,] -0.018892737 0.4176082 -0.9341080 0.9094812 0.532
## [7,] -0.055467763 0.4250969 -1.0942571 0.8192562 0.490
## [8,] -0.007563029 0.4117104 -0.9800507 0.8741712 0.548
## [9,] -0.032509036 0.4138104 -0.9341025 0.9188171 0.498
## [10,] -0.050876901 0.4142152 -1.0858341 0.8606469 0.548
## [11,] -0.043178598 0.4017999 -1.0177024 0.8189672 0.540
## [12,] -0.049915027 0.4154000 -0.9444616 0.8452286 0.518
## [13,] -0.044566791 0.4348905 -0.9736428 0.9201996 0.542
## [14,] -0.053132877 0.4085925 -0.9758844 0.8820159 0.482
## [15,] -0.046278234 0.3940232 -0.9180320 0.8092658 0.500
## [16,] -0.049514433 0.4578257 -1.1040893 0.9481966 0.534
## [17,] -0.070271101 0.4259853 -1.0340079 0.8111130 0.532
## [18,] -0.042616463 0.3950922 -1.0527976 0.8091935 0.482
## [19,] -0.039944456 0.4197640 -1.0431318 0.8395898 0.502plot(fit.npb.29$beta[,1], type = "l")
plot(fit.npb.29$beta[,2], type = "l")
plot(fit.npb.29$beta[,13], type = "l")
Same priors, but with scaleY = F
priors.npb.29a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 10e5, sig2inv.mu1 = 10e5)
fit.npb.29a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = FALSE,
priors = priors.npb.29a, interact = F)
npb.sum.29a <- summary(fit.npb.29a)
npb.sum.29a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.0000029929547 0.0008091636 -0.001938430 0.001936734 0.522
## [2,] 0.0000348936641 0.0007978201 -0.001778601 0.001874876 0.534
## [3,] 0.0000050700913 0.0008528772 -0.001883928 0.001829383 0.502
## [4,] -0.0000006592784 0.0008320206 -0.001783612 0.002043350 0.496
## [5,] 0.0000070091150 0.0007874591 -0.001657262 0.001945618 0.496
## [6,] -0.0000142948030 0.0008108816 -0.001693869 0.001954015 0.512
## [7,] 0.0000299715926 0.0008253950 -0.001715913 0.002258297 0.506
## [8,] 0.0000180048470 0.0008612665 -0.002026700 0.001880374 0.524
## [9,] 0.0000777242366 0.0007818032 -0.001610301 0.002026431 0.480
## [10,] -0.0000037155238 0.0008344709 -0.001948587 0.001890875 0.544
## [11,] 0.0000539999709 0.0008191653 -0.001793009 0.001999934 0.506
## [12,] -0.0000143277900 0.0007643464 -0.001657805 0.001806774 0.476
## [13,] 0.0000019952710 0.0008151975 -0.001850042 0.001889987 0.486
## [14,] -0.0000099848829 0.0008049083 -0.001912550 0.001996068 0.496
## [15,] 0.0000398545360 0.0007695658 -0.001654104 0.001944979 0.506
## [16,] 0.0000044317518 0.0008533496 -0.001914691 0.001868051 0.542
## [17,] 0.0000318213211 0.0007992530 -0.001655477 0.001985471 0.492
## [18,] 0.0000436334011 0.0007482470 -0.001591704 0.001803765 0.468
## [19,] 0.0000257047232 0.0007601775 -0.001762099 0.001869755 0.470plot(fit.npb.29a$beta[,1], type = "l")
plot(fit.npb.29a$beta[,2], type = "l")
plot(fit.npb.29a$beta[,13], type = "l")
4.1.7 Set alpha.pi and beta.pi to 5, set a.phi1 = 1, and adjust sig2inv.mu1
4.1.7.1 Try making a.phi1 = 1 and sig2inv.mu1 = 10
priors.npb.30 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 10)
fit.npb.30 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.30, interact = F)
npb.sum.30 <- summary(fit.npb.30)
npb.sum.30$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.7228556 6.185898 -20.101281 1.8610950 0.200
## [2,] -19.0050772 20.190536 -65.928054 0.0000000 0.684
## [3,] -1.3136830 5.208024 -19.738044 1.6554852 0.172
## [4,] -3.1161119 8.173765 -28.954626 0.2456618 0.236
## [5,] -0.4698456 4.195973 -11.667282 3.9323218 0.144
## [6,] -0.4558571 5.069872 -12.863175 6.0642046 0.160
## [7,] 0.3000572 6.026863 -8.381134 17.8194398 0.154
## [8,] 2.1301516 10.686663 -8.617136 41.4826373 0.194
## [9,] 0.2229029 6.518747 -9.693930 10.3608104 0.144
## [10,] -1.3429782 5.071480 -16.423261 1.5875183 0.178
## [11,] -3.9985155 9.784726 -35.178073 2.2364369 0.278
## [12,] -2.7480595 7.784423 -27.773997 1.6691937 0.242
## [13,] -4.3861962 9.853896 -35.243232 0.0000000 0.288
## [14,] -0.6861690 6.209898 -14.500452 2.4538156 0.170
## [15,] -3.5042988 8.989247 -32.428574 0.2456618 0.258
## [16,] -3.9430599 9.088170 -32.402623 0.4364748 0.284
## [17,] -9.4313789 16.160430 -50.994577 0.4967513 0.420
## [18,] -2.0118603 6.588875 -19.891292 0.0000000 0.212
## [19,] -1.6924124 6.275389 -19.649915 0.8215062 0.170plot(fit.npb.30$beta[,1], type = "l")
plot(fit.npb.30$beta[,2], type = "l")
plot(fit.npb.30$beta[,13], type = "l")
Same priors, but with scaleY = F
priors.npb.30a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 10)
fit.npb.30a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = FALSE,
priors = priors.npb.30a, interact = F)
npb.sum.30a <- summary(fit.npb.30a)
npb.sum.30a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] 0.074119118 1.363934 -2.308498 2.756957 0.468
## [2,] 0.010494643 1.326295 -2.420825 2.808571 0.512
## [3,] 0.063236533 1.262993 -2.227545 2.588655 0.482
## [4,] 0.004503723 1.485491 -2.346604 2.475782 0.488
## [5,] -0.018773342 1.337000 -2.403715 2.326143 0.510
## [6,] -0.014432961 1.298491 -2.198557 2.679880 0.476
## [7,] 0.098224373 1.286374 -1.921173 2.803585 0.476
## [8,] -0.046045737 1.194057 -2.256074 2.550880 0.502
## [9,] 0.060598256 1.322770 -2.188775 2.812985 0.516
## [10,] 0.103873445 1.319422 -2.204537 2.905311 0.488
## [11,] 0.067135740 1.344225 -2.214774 1.920478 0.530
## [12,] 0.064489504 1.228626 -2.179235 2.507171 0.496
## [13,] -0.038101943 1.473007 -3.643061 2.745544 0.480
## [14,] 0.056993085 1.264079 -2.323791 3.011037 0.496
## [15,] 0.043757160 1.258120 -2.222910 2.656063 0.504
## [16,] 0.121804354 1.072164 -1.878659 2.567970 0.468
## [17,] -0.020926185 1.240334 -3.079505 2.507171 0.492
## [18,] 0.072878661 1.423697 -2.532435 2.588655 0.480
## [19,] -0.032569195 1.034985 -1.843805 1.854254 0.456plot(fit.npb.30a$beta[,1], type = "l")
plot(fit.npb.30a$beta[,2], type = "l")
plot(fit.npb.30a$beta[,13], type = "l")
4.1.7.2 Try making a.phi1 = 1 and sig2inv.mu1 = 100
priors.npb.31 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 100)
fit.npb.31 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.31, interact = F)
npb.sum.31 <- summary(fit.npb.31)
npb.sum.31$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -2.30889214 6.784902 -24.771436 0.000000000 0.200
## [2,] -20.92802151 20.179046 -62.805035 0.000000000 0.732
## [3,] -1.08965652 4.397055 -12.799255 0.000000000 0.158
## [4,] -3.37844361 8.694089 -32.124019 0.003170527 0.252
## [5,] -0.93522724 4.686011 -13.243451 0.188809267 0.156
## [6,] -0.54554790 5.140308 -12.597600 6.688217638 0.164
## [7,] -0.19526708 5.578431 -11.915487 12.222610797 0.184
## [8,] 1.75217377 9.034742 -7.926483 33.065635224 0.174
## [9,] 0.09622656 5.085413 -8.919991 12.543905370 0.146
## [10,] -1.07175896 6.258691 -16.323661 6.033725382 0.192
## [11,] -4.87259495 10.790576 -36.999806 0.003170527 0.310
## [12,] -2.65268645 7.873357 -25.620569 0.000000000 0.238
## [13,] -4.77325818 10.868420 -37.171170 0.038797733 0.314
## [14,] -0.70896380 5.283107 -12.090274 5.821940112 0.172
## [15,] -3.02747698 8.213266 -28.371188 0.538732726 0.266
## [16,] -3.11612392 8.292508 -28.674739 0.000000000 0.260
## [17,] -11.12370533 17.011735 -56.854870 0.000000000 0.454
## [18,] -2.24537034 7.309293 -26.808510 0.108378804 0.214
## [19,] -2.18288162 6.462647 -24.682852 0.038797733 0.220plot(fit.npb.31$beta[,1], type = "l")
plot(fit.npb.31$beta[,2], type = "l")
plot(fit.npb.31$beta[,13], type = "l")
Same priors, but with scaleY = F
priors.npb.31a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 100)
fit.npb.31a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = FALSE,
priors = priors.npb.31a, interact = F)
npb.sum.31a <- summary(fit.npb.31a)
npb.sum.31a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -0.11200176 1.845975 -3.594599 2.933795 0.524
## [2,] -0.06084745 1.446065 -3.230574 2.402961 0.544
## [3,] -0.03186257 1.609605 -3.502165 2.666465 0.494
## [4,] -0.05796699 1.565341 -3.072019 2.393704 0.542
## [5,] -0.03733648 1.353984 -2.690057 2.838679 0.500
## [6,] -0.10830473 1.893346 -5.062828 2.871385 0.532
## [7,] -0.09578467 1.428884 -3.525617 2.485758 0.504
## [8,] -0.10651158 1.903310 -3.568060 2.641535 0.518
## [9,] -0.07269353 1.475613 -3.204855 2.616656 0.474
## [10,] -0.16545997 1.853086 -3.372041 2.382987 0.518
## [11,] -0.18876057 1.724617 -3.552851 2.094670 0.478
## [12,] -0.09518493 1.543222 -3.608179 2.813356 0.550
## [13,] -0.08358490 1.709686 -3.763812 2.924885 0.526
## [14,] 0.01514462 1.481139 -3.398926 3.719411 0.490
## [15,] -0.16326693 1.731134 -3.794918 2.855899 0.508
## [16,] -0.02413426 1.722385 -3.230866 3.497753 0.518
## [17,] -0.10137685 1.957256 -3.491128 2.648940 0.482
## [18,] 0.02946311 1.412729 -2.854176 2.989497 0.522
## [19,] -0.02930316 1.780401 -3.267473 2.635170 0.524plot(fit.npb.31a$beta[,1], type = "l")
plot(fit.npb.31a$beta[,2], type = "l")
plot(fit.npb.31a$beta[,13], type = "l")
4.1.7.3 Try making a.phi1 = 1 and sig2inv.mu1 = 10e5
priors.npb.32 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 10e5)
fit.npb.32 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb.32, interact = F)
npb.sum.32 <- summary(fit.npb.32)
npb.sum.32$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.38446846 7.022531 -20.332156 10.17422656 0.226
## [2,] -20.62264930 20.077258 -65.089736 0.08685412 0.720
## [3,] -1.06435183 4.917391 -14.901425 2.53796662 0.166
## [4,] -3.06679682 8.528214 -27.692337 2.01332470 0.250
## [5,] -0.84960024 5.017413 -16.335625 5.41330050 0.146
## [6,] -0.38439458 5.549323 -13.393940 8.34579300 0.164
## [7,] 0.03984075 5.288480 -12.413025 10.51676007 0.148
## [8,] 2.37295219 10.314854 -6.581860 39.89358420 0.170
## [9,] 0.54030411 6.489239 -9.616982 16.48935580 0.166
## [10,] -0.64553010 5.324006 -14.722058 7.73864664 0.168
## [11,] -5.08187604 11.058515 -36.058333 2.46512327 0.308
## [12,] -3.02452916 8.638986 -31.762681 0.98918597 0.226
## [13,] -4.79160611 10.602028 -37.257283 0.42522579 0.290
## [14,] -0.47361111 5.411474 -12.858097 11.24888674 0.170
## [15,] -3.14471481 8.448053 -30.112915 2.14570556 0.232
## [16,] -2.96886714 9.279894 -34.730379 4.91918466 0.250
## [17,] -10.16885803 16.029957 -51.528538 0.26003816 0.430
## [18,] -2.44958126 8.455490 -28.831080 4.40927059 0.208
## [19,] -2.15767429 7.679553 -24.673937 7.73864664 0.214plot(fit.npb.32$beta[,1], type = "l")
plot(fit.npb.32$beta[,2], type = "l")
plot(fit.npb.32$beta[,13], type = "l")
Same priors, but with scaleY = F
priors.npb.32a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
a.phi1 = 1, sig2inv.mu1 = 10e5)
fit.npb.32a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = FALSE,
priors = priors.npb.32a, interact = F)
npb.sum.32a <- summary(fit.npb.32a)
npb.sum.32a$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] 0.025209362 1.294861 -2.347621 2.420512 0.514
## [2,] -0.049016935 1.341006 -2.409473 2.263052 0.494
## [3,] 0.028896975 1.381161 -1.960469 2.326273 0.518
## [4,] -0.038418735 1.219774 -2.483217 1.975696 0.518
## [5,] -0.023985814 1.355683 -2.392673 2.356443 0.474
## [6,] 0.030791522 1.302866 -2.599089 2.388541 0.478
## [7,] -0.152837433 1.720482 -2.534500 1.739969 0.486
## [8,] -0.040706452 1.508916 -2.671380 2.258912 0.492
## [9,] 0.092356015 1.371136 -2.119311 2.417871 0.480
## [10,] -0.075665075 1.216456 -2.298597 2.012301 0.496
## [11,] 0.029002313 1.263324 -2.564436 2.791453 0.486
## [12,] -0.099399846 1.027492 -2.148132 2.105168 0.436
## [13,] 0.152498936 1.428497 -2.158852 3.180683 0.498
## [14,] 0.062607859 1.158282 -1.939363 2.221776 0.482
## [15,] -0.054083611 1.694672 -2.491200 2.768224 0.496
## [16,] -0.064899159 1.629433 -3.019469 2.619183 0.514
## [17,] -0.099680901 1.190646 -2.386633 1.873959 0.502
## [18,] -0.008919667 1.258053 -2.299866 2.293204 0.496
## [19,] 0.092523129 1.529017 -2.054268 2.607070 0.494plot(fit.npb.32a$beta[,1], type = "l")
plot(fit.npb.32a$beta[,2], type = "l")
plot(fit.npb.32a$beta[,13], type = "l")
4.2 Fit the NPB model
Below I’ve used the 24th set of priors and set scaleY = T
The priors are as follows: r priors.npb.24
Note that this version of the model does not include gest_age_w.
priors.npb <- priors.npb.24
fit.npb <- npb(niter = 10000, nburn = 2500, X = X.scaled, Y = Y, W = W.scaled2,
scaleY = TRUE,
priors = priors.npb, interact = TRUE, XWinteract = TRUE)
save(fit.npb, file = here::here("Results", "NPB_Birth_Weight.rdata"))
# load(here::here("Results", "NPB_Birth_Weight.rdata"))
npb.sum <- summary(fit.npb)4.2.1 First, main effect regression coefficients with PIPs
npb.sum$main.effects## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] -1.8661265 6.914319 -22.957048 6.049375 0.2281333
## [2,] -24.9299099 21.166516 -69.288505 0.000000 0.7853333
## [3,] -1.1551286 6.033682 -18.834106 8.098385 0.2098667
## [4,] -3.1530965 9.092592 -32.056663 3.639703 0.2676000
## [5,] -0.5165156 5.370372 -14.232302 9.672763 0.1869333
## [6,] -0.2570868 5.721849 -14.272370 12.276842 0.1908000
## [7,] 0.5523885 7.135228 -11.182066 19.879487 0.1908000
## [8,] 3.8083493 12.832219 -7.874556 47.771370 0.2461333
## [9,] 0.7159346 7.048218 -10.640137 20.509450 0.1908000
## [10,] -0.8666888 5.859006 -17.095504 7.869133 0.1969333
## [11,] -4.7334937 11.084603 -39.038464 2.605766 0.3028000
## [12,] -2.3870316 8.563241 -29.073269 6.775860 0.2480000
## [13,] -5.3488304 11.913380 -42.129727 1.969438 0.3268000
## [14,] -0.4622358 6.634494 -15.841082 12.762775 0.2030667
## [15,] -3.3172436 8.902109 -31.384447 3.383887 0.2705333
## [16,] -3.0570276 9.676074 -32.593852 7.114270 0.2761333
## [17,] -12.3144543 18.075048 -57.590304 0.000000 0.4890667
## [18,] -1.7492751 7.631214 -23.981519 7.848185 0.2302667
## [19,] -2.4616124 9.410943 -31.212532 7.106345 0.2552000selected_exp <- which(npb.sum$main.effects[,5] >= 0.5)
selected_exp## [1] 2#' Which variables are these?
exp_names <- colnames(X)[selected_exp]
exp_names## [1] "mean_o3"4.2.2 Interactions
Next, all of the interactions between exposures or between exposures and covariates
npb.sum$interactions## Posterior Mean SD 95% CI Lower 95% CI Upper PIP
## [1,] 0.00207935185 0.3922357 0 0 0.0026666667
## [2,] -0.01092256325 0.8419317 0 0 0.0030666667
## [3,] 0.00406659348 0.2943407 0 0 0.0012000000
## [4,] 0.02102702601 0.8433957 0 0 0.0018666667
## [5,] -0.00120801985 0.2514967 0 0 0.0012000000
## [6,] -0.00434257270 0.5189284 0 0 0.0025333333
## [7,] 0.00162709402 0.2961235 0 0 0.0013333333
## [8,] -0.00828639526 0.7148152 0 0 0.0025333333
## [9,] -0.04316211873 1.1765333 0 0 0.0040000000
## [10,] -0.01240167260 0.4834403 0 0 0.0017333333
## [11,] -0.00816117136 0.6874690 0 0 0.0024000000
## [12,] -0.00174916587 0.3175548 0 0 0.0020000000
## [13,] -0.03365644311 0.9396474 0 0 0.0025333333
## [14,] -0.00081702879 0.3060190 0 0 0.0014666667
## [15,] -0.01647023314 0.7112123 0 0 0.0021333333
## [16,] -0.00026800254 0.2458204 0 0 0.0014666667
## [17,] -0.01152744851 0.3727665 0 0 0.0025333333
## [18,] -0.00362124603 0.4932116 0 0 0.0025333333
## [19,] -0.02240938075 0.8128801 0 0 0.0028000000
## [20,] 0.01113662913 0.4821056 0 0 0.0016000000
## [21,] -0.02836486280 0.8092080 0 0 0.0025333333
## [22,] -0.10368334545 1.9930455 0 0 0.0058666667
## [23,] -0.06378075968 1.4115292 0 0 0.0036000000
## [24,] -0.17107318957 2.4115089 0 0 0.0077333333
## [25,] -0.00986094243 0.5076520 0 0 0.0025333333
## [26,] -0.14250513955 2.2457617 0 0 0.0062666667
## [27,] -0.04747840190 1.2283511 0 0 0.0034666667
## [28,] 0.00978619120 0.4129487 0 0 0.0010666667
## [29,] -0.00351463226 0.3097896 0 0 0.0020000000
## [30,] -0.00738293457 0.3842999 0 0 0.0012000000
## [31,] 0.00486774189 0.3790955 0 0 0.0021333333
## [32,] 0.03384653345 1.0808633 0 0 0.0028000000
## [33,] 0.00654153167 0.3261972 0 0 0.0009333333
## [34,] 0.00625761431 0.4010372 0 0 0.0017333333
## [35,] 0.06025408379 1.5010229 0 0 0.0030666667
## [36,] 0.01155648563 0.4596850 0 0 0.0020000000
## [37,] 0.00177964328 0.4078900 0 0 0.0021333333
## [38,] -0.00193298294 0.3107592 0 0 0.0009333333
## [39,] -0.00573018811 0.4239365 0 0 0.0022666667
## [40,] 0.00082028627 0.5543440 0 0 0.0022666667
## [41,] 0.01204786388 0.4962271 0 0 0.0018666667
## [42,] -0.02178836559 0.9648496 0 0 0.0029333333
## [43,] -0.01684532461 0.9548668 0 0 0.0028000000
## [44,] -0.00675761965 0.3395086 0 0 0.0028000000
## [45,] -0.00613913712 0.4891281 0 0 0.0016000000
## [46,] -0.00762962299 0.4389237 0 0 0.0016000000
## [47,] -0.00663410436 0.3756887 0 0 0.0016000000
## [48,] 0.00669872724 0.3964017 0 0 0.0017333333
## [49,] -0.00034435030 0.1390577 0 0 0.0010666667
## [50,] 0.01833404806 0.8046079 0 0 0.0029333333
## [51,] 0.04158558967 1.4713239 0 0 0.0022666667
## [52,] 0.01122000351 0.5265119 0 0 0.0014666667
## [53,] 0.00617330700 0.3277792 0 0 0.0016000000
## [54,] -0.00319363424 0.2111496 0 0 0.0014666667
## [55,] 0.01188933563 0.5087609 0 0 0.0025333333
## [56,] 0.00167179752 0.3188443 0 0 0.0018666667
## [57,] -0.00953491420 0.3848878 0 0 0.0014666667
## [58,] 0.00932028287 0.6625620 0 0 0.0021333333
## [59,] -0.01507875291 0.5634662 0 0 0.0018666667
## [60,] -0.00903040186 0.3625436 0 0 0.0021333333
## [61,] 0.00009563128 0.1890493 0 0 0.0012000000
## [62,] -0.00106960931 0.2316240 0 0 0.0014666667
## [63,] -0.00010282415 0.4312542 0 0 0.0024000000
## [64,] -0.00819730528 0.4278985 0 0 0.0017333333
## [65,] -0.00271505765 0.5355645 0 0 0.0024000000
## [66,] 0.00696295275 0.3505754 0 0 0.0025333333
## [67,] 0.00886277979 0.6603891 0 0 0.0021333333
## [68,] 0.01227569826 0.4671397 0 0 0.0021333333
## [69,] 0.02017797812 0.6489430 0 0 0.0033333333
## [70,] 0.03460275885 1.0235756 0 0 0.0025333333
## [71,] -0.00531811777 0.3803452 0 0 0.0017333333
## [72,] -0.02482249166 0.7397177 0 0 0.0025333333
## [73,] 0.00136461852 0.3020055 0 0 0.0021333333
## [74,] -0.01021762302 0.5697642 0 0 0.0025333333
## [75,] -0.02306432047 0.7495557 0 0 0.0030666667
## [76,] -0.02784400404 0.8632075 0 0 0.0030666667
## [77,] 0.00790505529 0.5734143 0 0 0.0025333333
## [78,] -0.00056101225 0.2453785 0 0 0.0021333333
## [79,] 0.01096704290 0.4862492 0 0 0.0034666667
## [80,] 0.00635613678 0.4227833 0 0 0.0012000000
## [81,] -0.00040256004 0.2218477 0 0 0.0014666667
## [82,] 0.00515834614 0.5372605 0 0 0.0020000000
## [83,] -0.00007035086 0.4211832 0 0 0.0022666667
## [84,] 0.00610592089 0.4112330 0 0 0.0026666667
## [85,] -0.00296153989 0.4252058 0 0 0.0025333333
## [86,] 0.00233603166 0.4785649 0 0 0.0018666667
## [87,] 0.00950760536 0.4375328 0 0 0.0016000000
## [88,] -0.00319741519 0.2553187 0 0 0.0012000000
## [89,] 0.00279567956 0.2591166 0 0 0.0012000000
## [90,] -0.01095113162 0.7119509 0 0 0.0018666667
## [91,] 0.01754758253 0.7789657 0 0 0.0021333333
## [92,] -0.00248372582 0.3042499 0 0 0.0018666667
## [93,] 0.01056954134 0.6007759 0 0 0.0026666667
## [94,] 0.01354972302 0.8107644 0 0 0.0020000000
## [95,] 0.00610162425 0.2657757 0 0 0.0018666667
## [96,] 0.00550132255 0.3047765 0 0 0.0008000000
## [97,] -0.00328594477 0.3239822 0 0 0.0013333333
## [98,] 0.00824007988 0.5886503 0 0 0.0025333333
## [99,] 0.00377064657 0.3345422 0 0 0.0017333333
## [100,] 0.00268087629 0.4923001 0 0 0.0017333333
## [101,] 0.00361759755 0.2664276 0 0 0.0016000000
## [102,] -0.00125785535 0.1884398 0 0 0.0020000000
## [103,] 0.00277632737 0.4470401 0 0 0.0022666667
## [104,] -0.01170916706 0.7070968 0 0 0.0024000000
## [105,] -0.00950454374 0.5092766 0 0 0.0014666667
## [106,] 0.00991839707 0.4102977 0 0 0.0018666667
## [107,] 0.00760498155 0.7110995 0 0 0.0030666667
## [108,] -0.00717883556 0.5464990 0 0 0.0029333333
## [109,] -0.00824565656 0.6836161 0 0 0.0024000000
## [110,] 0.00751171461 0.3115849 0 0 0.0013333333
## [111,] -0.00055621021 0.2717771 0 0 0.0021333333
## [112,] 0.00084782172 0.3858959 0 0 0.0018666667
## [113,] -0.01374735320 0.5549056 0 0 0.0024000000
## [114,] -0.00030905418 0.4253479 0 0 0.0016000000
## [115,] 0.00044446584 0.3504592 0 0 0.0020000000
## [116,] 0.00892639020 0.3272707 0 0 0.0020000000
## [117,] 0.03063198668 0.8941522 0 0 0.0024000000
## [118,] 0.05853421103 1.3669008 0 0 0.0036000000
## [119,] -0.01799978263 0.5900237 0 0 0.0024000000
## [120,] -0.01381650457 0.6588336 0 0 0.0024000000
## [121,] -0.00801616383 0.4098493 0 0 0.0012000000
## [122,] 0.00231720140 0.2151038 0 0 0.0014666667
## [123,] -0.02841886236 0.8972683 0 0 0.0034666667
## [124,] -0.02196784884 0.8551798 0 0 0.0024000000
## [125,] -0.01820907688 0.6794758 0 0 0.0021333333
## [126,] 0.00191917596 0.4796340 0 0 0.0028000000
## [127,] -0.05492967323 1.1487590 0 0 0.0037333333
## [128,] 0.00075666617 0.4204341 0 0 0.0024000000
## [129,] -0.00192662113 0.3921484 0 0 0.0017333333
## [130,] -0.00894344962 0.5742215 0 0 0.0017333333
## [131,] -0.01273999395 0.5222802 0 0 0.0012000000
## [132,] -0.00101176390 0.6148757 0 0 0.0016000000
## [133,] -0.01271177986 0.7720026 0 0 0.0018666667
## [134,] 0.00198568213 0.8165910 0 0 0.0024000000
## [135,] -0.00614984337 0.5660170 0 0 0.0020000000
## [136,] 0.00977359504 0.6418780 0 0 0.0020000000
## [137,] 0.01685598053 0.5467248 0 0 0.0020000000
## [138,] 0.00478643032 0.6623801 0 0 0.0022666667
## [139,] -0.00148768633 0.1901172 0 0 0.0008000000
## [140,] 0.03684251473 1.0686060 0 0 0.0025333333
## [141,] 0.02237133934 0.7820852 0 0 0.0025333333
## [142,] 0.08334718299 2.2104629 0 0 0.0041333333
## [143,] 0.01670596339 0.7741902 0 0 0.0026666667
## [144,] -0.00646702378 0.3409709 0 0 0.0010666667
## [145,] -0.01787038539 0.6178606 0 0 0.0020000000
## [146,] -0.01210710502 0.5491503 0 0 0.0017333333
## [147,] -0.01766597014 0.7662021 0 0 0.0021333333
## [148,] -0.13610005814 2.1474472 0 0 0.0064000000
## [149,] -0.03446115081 0.9651905 0 0 0.0024000000
## [150,] -0.05639773821 1.2668824 0 0 0.0033333333
## [151,] -0.02135982801 0.6118892 0 0 0.0017333333
## [152,] -0.00711709735 0.2611309 0 0 0.0016000000
## [153,] -0.00271554287 0.1810662 0 0 0.0006666667
## [154,] -0.01551848705 0.6425994 0 0 0.0029333333
## [155,] 0.00231839316 0.4860063 0 0 0.0024000000
## [156,] -0.05187495157 1.1410308 0 0 0.0034666667
## [157,] -0.00092254702 0.1006495 0 0 0.0014666667
## [158,] -0.00349315349 0.2194937 0 0 0.0017333333
## [159,] 0.00638799743 0.3075779 0 0 0.0016000000
## [160,] -0.00319303611 0.5067947 0 0 0.0017333333
## [161,] 0.00187973045 0.2029551 0 0 0.0010666667
## [162,] -0.00083055311 0.4175919 0 0 0.0014666667
## [163,] 0.00722467650 0.4821929 0 0 0.0022666667
## [164,] 0.01699282389 0.7166051 0 0 0.0028000000
## [165,] -0.00236839311 0.4269201 0 0 0.0028000000
## [166,] -0.00704964679 0.2913808 0 0 0.0010666667
## [167,] 0.00231751600 0.2206739 0 0 0.0016000000
## [168,] 0.00163668294 0.2172364 0 0 0.0010666667
## [169,] -0.01664527911 0.6274140 0 0 0.0020000000
## [170,] -0.00899307613 0.3302830 0 0 0.0018666667
## [171,] 0.00280111222 0.2434321 0 0 0.0014666667
## [172,] -0.00081297778 0.6886768 0 0 0.0032000000
## [173,] -0.03431340239 1.1892105 0 0 0.0032000000
## [174,] 0.01020203271 1.3276948 0 0 0.0034666667
## [175,] -0.00306244740 1.1023503 0 0 0.0025333333
## [176,] -0.04166310150 1.5189742 0 0 0.0025333333
## [177,] -0.04865627529 1.3699596 0 0 0.0029333333
## [178,] 0.02272181788 1.1637376 0 0 0.0021333333
## [179,] -0.06124256069 1.8908378 0 0 0.0033333333
## [180,] -0.09466697984 2.0758565 0 0 0.0037333333
## [181,] 0.00520553456 1.0137846 0 0 0.0026666667
## [182,] 0.00038624600 0.2136883 0 0 0.0016000000
## [183,] -0.00775767903 1.4737020 0 0 0.0028000000
## [184,] -0.00737470658 0.6535756 0 0 0.0024000000
## [185,] 0.01983007156 0.5953587 0 0 0.0028000000
## [186,] -0.00808493610 0.5810289 0 0 0.0028000000
## [187,] -0.01034594011 0.6195578 0 0 0.0017333333
## [188,] -0.01815461465 0.9762409 0 0 0.0020000000
## [189,] -0.02674016229 1.1013299 0 0 0.0020000000
## [190,] -0.01556922389 1.2275186 0 0 0.0024000000
## [191,] -0.09440391118 2.2734959 0 0 0.0037333333
## [192,] -0.11541966605 2.8579405 0 0 0.0042666667
## [193,] -0.06505030255 1.8299151 0 0 0.0030666667
## [194,] -0.01912226918 0.9089227 0 0 0.0026666667
## [195,] 0.02722410942 2.6070112 0 0 0.0030666667
## [196,] -0.06220270214 1.9149677 0 0 0.0028000000
## [197,] -0.05534092807 2.1367940 0 0 0.0030666667
## [198,] -0.00754190509 0.4270707 0 0 0.0016000000
## [199,] -0.03940616146 1.7275463 0 0 0.0029333333
## [200,] -0.38716295306 5.5544717 0 0 0.0082666667
## [201,] 0.04122867996 1.2606507 0 0 0.0028000000
## [202,] -0.01060857689 0.4979391 0 0 0.0021333333
## [203,] 0.00162961556 0.4884729 0 0 0.0021333333
## [204,] 0.01036658606 1.1853954 0 0 0.0025333333
## [205,] -0.00338327579 1.1089221 0 0 0.0020000000
## [206,] 0.00435476674 1.3530703 0 0 0.0028000000
## [207,] 0.05356144388 2.6682333 0 0 0.0026666667
## [208,] -0.05672995038 2.5429039 0 0 0.0028000000
## [209,] 0.04243628741 1.7509027 0 0 0.0024000000
## [210,] -0.02046091494 0.7706384 0 0 0.0028000000
## [211,] -0.02405431603 1.7130633 0 0 0.0032000000
## [212,] -0.02422102191 0.9172878 0 0 0.0030666667
## [213,] 0.01101646059 1.4250734 0 0 0.0020000000
## [214,] -0.00063976827 0.3329916 0 0 0.0018666667
## [215,] -0.00043308647 0.7750810 0 0 0.0013333333
## [216,] -0.00164320983 0.4185761 0 0 0.0022666667
## [217,] 0.00744330952 0.4244182 0 0 0.0014666667
## [218,] 0.00155443218 0.4258177 0 0 0.0013333333
## [219,] -0.00215161030 0.6208105 0 0 0.0018666667
## [220,] -0.02766362758 0.9629175 0 0 0.0029333333
## [221,] -0.00556598363 0.7027008 0 0 0.0032000000
## [222,] -0.02051517955 1.1680229 0 0 0.0020000000
## [223,] 0.02136577403 2.0597233 0 0 0.0021333333
## [224,] 0.01420880434 0.9716482 0 0 0.0026666667
## [225,] -0.03947772384 1.0719718 0 0 0.0029333333
## [226,] -0.00990273047 0.7471768 0 0 0.0018666667
## [227,] -0.02368929216 2.1886330 0 0 0.0025333333
## [228,] -0.02507830105 0.8714135 0 0 0.0029333333
## [229,] 0.00191519565 0.4649630 0 0 0.0021333333
## [230,] -0.02356780432 0.7582051 0 0 0.0025333333
## [231,] 0.01969205293 1.0462867 0 0 0.0028000000
## [232,] 0.01062621297 0.6339565 0 0 0.0017333333
## [233,] -0.04767319049 1.1191493 0 0 0.0034666667
## [234,] -0.00624219615 0.5686079 0 0 0.0013333333
## [235,] -0.07297478542 1.7677986 0 0 0.0037333333
## [236,] 0.09308434432 3.1620765 0 0 0.0037333333
## [237,] 0.01367961526 0.8796297 0 0 0.0017333333
## [238,] -0.03100237972 1.1793412 0 0 0.0032000000
## [239,] 0.18407676227 4.1054085 0 0 0.0045333333
## [240,] 0.01925506764 1.4312522 0 0 0.0020000000
## [241,] -0.00332800596 0.7648545 0 0 0.0026666667
## [242,] -0.02034019953 1.0077404 0 0 0.0034666667
## [243,] -0.01151928020 0.9615075 0 0 0.0026666667
## [244,] -0.08346867671 1.8703253 0 0 0.0048000000
## [245,] 0.01068482020 0.8011117 0 0 0.0014666667
## [246,] -0.12704828267 2.0904095 0 0 0.0058666667
## [247,] -0.01582609695 0.8666925 0 0 0.0022666667
## [248,] -0.00664341987 0.5205422 0 0 0.0025333333
## [249,] -0.00768390617 0.3883377 0 0 0.0018666667
## [250,] 0.01747066456 0.6650394 0 0 0.0018666667
## [251,] -0.00964438519 0.5120514 0 0 0.0026666667
## [252,] 0.01272525738 1.7705447 0 0 0.0016000000
## [253,] -0.01949159208 0.7819220 0 0 0.0025333333
## [254,] -0.01046251581 1.0086585 0 0 0.0014666667
## [255,] 0.23065547591 6.2896192 0 0 0.0032000000
## [256,] -0.06016821592 2.2611802 0 0 0.0036000000
## [257,] 0.01317283119 0.7031169 0 0 0.0026666667
## [258,] -0.00460046878 0.4615307 0 0 0.0025333333
## [259,] -0.10953812678 4.3456181 0 0 0.0034666667
## [260,] -0.00779575824 0.5559891 0 0 0.0020000000
## [261,] 0.00383455567 1.2674115 0 0 0.0022666667
## [262,] -0.00273202064 0.2961800 0 0 0.0020000000
## [263,] -0.02486766639 1.5239684 0 0 0.0029333333
## [264,] -0.03221595933 1.4038361 0 0 0.0021333333
## [265,] 0.00584595257 0.6471392 0 0 0.0022666667
## [266,] -0.01288032877 0.5282785 0 0 0.0022666667
## [267,] 0.00958446685 0.4937792 0 0 0.0010666667
## [268,] -0.00051543682 0.8231363 0 0 0.0017333333
## [269,] 0.00296068457 0.4619030 0 0 0.0022666667
## [270,] -0.01110336188 1.0910313 0 0 0.0021333333
## [271,] 0.10569172239 3.2441184 0 0 0.0037333333
## [272,] -0.00500407114 0.7799976 0 0 0.0013333333
## [273,] 0.01348582480 1.0429696 0 0 0.0026666667
## [274,] 0.00214440035 1.0787527 0 0 0.0026666667
## [275,] -0.03961722952 1.8730041 0 0 0.0033333333
## [276,] 0.01139647608 0.8074230 0 0 0.0024000000
## [277,] -0.01103854317 0.6308667 0 0 0.0017333333
## [278,] -0.01205514139 0.6580385 0 0 0.0029333333
## [279,] 0.00257183737 1.2283511 0 0 0.0029333333
## [280,] -0.09736940558 2.4403745 0 0 0.0038666667
## [281,] 0.01163660297 0.4175403 0 0 0.0020000000
## [282,] -0.06649922248 1.5648285 0 0 0.0037333333
## [283,] 0.01098834837 0.8059703 0 0 0.0028000000
## [284,] -0.00052475651 0.6636149 0 0 0.0024000000
## [285,] -0.00718199968 0.9689494 0 0 0.0022666667
## [286,] 0.05676040446 2.5439325 0 0 0.0030666667
## [287,] -0.00077446957 1.2032435 0 0 0.0022666667
## [288,] 0.03798122587 1.6691974 0 0 0.0018666667
## [289,] -0.01888351315 0.7508118 0 0 0.0025333333
## [290,] 0.06999739199 1.8038485 0 0 0.0038666667
## [291,] -0.03747363234 1.5030872 0 0 0.0025333333
## [292,] 0.05936526073 1.9413802 0 0 0.0029333333
## [293,] 0.01886938863 1.4535352 0 0 0.0026666667
## [294,] 0.00174599762 0.4802649 0 0 0.0021333333
## [295,] 0.02146633259 1.4784986 0 0 0.0017333333
## [296,] -0.01003938128 0.8981341 0 0 0.0021333333
## [297,] 0.00246918798 0.4344574 0 0 0.0018666667
## [298,] -0.00214915651 0.4119082 0 0 0.0017333333
## [299,] 0.03585327099 1.4999970 0 0 0.0021333333
## [300,] -0.01204974994 0.7223358 0 0 0.0033333333
## [301,] 0.10001951896 3.3922676 0 0 0.0038666667
## [302,] 0.00795935423 0.6471900 0 0 0.0020000000
## [303,] 0.00756847533 1.0036607 0 0 0.0030666667
## [304,] 0.00462243457 0.8571579 0 0 0.0032000000
## [305,] 0.00618317451 0.7161878 0 0 0.0025333333
## [306,] -0.00308465839 0.8810382 0 0 0.0020000000
## [307,] -0.03995431964 2.9813258 0 0 0.0021333333
## [308,] 0.03570797249 1.2577194 0 0 0.0032000000
## [309,] -0.00590692519 0.4147279 0 0 0.0020000000
## [310,] -0.00587217822 0.2587566 0 0 0.0017333333
## [311,] -0.00419693225 0.8878103 0 0 0.0021333333
## [312,] -0.00358781651 0.7196047 0 0 0.0020000000
## [313,] 0.00009117405 0.2979145 0 0 0.0021333333
## [314,] -0.00794283133 0.6021553 0 0 0.0021333333
## [315,] 0.01845327463 0.8951495 0 0 0.0022666667
## [316,] -0.02301533698 0.9855940 0 0 0.0016000000
## [317,] -0.02322966334 0.9100376 0 0 0.0028000000
## [318,] 0.00878496706 0.8142470 0 0 0.0022666667
## [319,] 0.07016490094 2.9730100 0 0 0.0033333333
## [320,] -0.02754680319 0.9148212 0 0 0.0029333333
## [321,] -0.00405799753 0.6138842 0 0 0.0030666667
## [322,] -0.01381030436 0.9571022 0 0 0.0024000000
## [323,] -0.35308212308 9.4437342 0 0 0.0042666667
## [324,] -0.00695767588 0.6764322 0 0 0.0017333333
## [325,] 0.00374329079 0.3925999 0 0 0.0021333333
## [326,] 0.00018037313 0.6355554 0 0 0.0020000000
## [327,] -0.02371853914 0.9696055 0 0 0.0025333333
## [328,] -0.06570402085 1.7301026 0 0 0.0033333333
## [329,] -0.00375723919 0.2468508 0 0 0.0017333333
## [330,] 0.00094764973 0.2418656 0 0 0.0013333333
## [331,] -0.00119904665 0.6752365 0 0 0.0018666667
## [332,] -0.01255210481 0.8980907 0 0 0.0017333333
## [333,] -0.02762448463 0.9563711 0 0 0.0022666667
## [334,] 0.01351343700 1.7914729 0 0 0.0026666667
## [335,] -0.02080827693 0.9893235 0 0 0.0026666667
## [336,] -0.07854615426 2.0733521 0 0 0.0028000000
## [337,] 0.01313104320 0.7623468 0 0 0.0024000000
## [338,] -0.00973182484 0.4518883 0 0 0.0013333333
## [339,] -0.00206754392 1.0372468 0 0 0.0030666667
## [340,] -0.01999001313 0.9746472 0 0 0.0029333333
## [341,] 0.00071573863 0.3526947 0 0 0.0018666667
## [342,] 0.00260855587 0.4130259 0 0 0.0022666667
## [343,] -0.01078974107 0.9436231 0 0 0.0024000000
## [344,] -0.09997933831 2.4714020 0 0 0.0036000000
## [345,] 0.00581780556 0.5477047 0 0 0.0020000000
## [346,] 0.03464659912 0.9728667 0 0 0.0032000000
## [347,] -0.02746845970 0.8049390 0 0 0.0030666667
## [348,] -0.05196948929 1.6847849 0 0 0.0036000000
## [349,] -0.00770447005 0.9568034 0 0 0.0024000000
## [350,] -0.03548639673 1.3588053 0 0 0.0030666667
## [351,] 0.01746270507 2.0771395 0 0 0.0018666667
## [352,] -0.00976209087 0.5721004 0 0 0.0018666667
## [353,] 0.01307056949 1.1395545 0 0 0.0021333333
## [354,] -0.07073266794 2.1072539 0 0 0.0042666667
## [355,] -0.00327447574 0.5864519 0 0 0.0028000000
## [356,] -0.00159596872 0.6706603 0 0 0.0020000000
## [357,] 0.00504552329 0.8321624 0 0 0.0024000000
## [358,] 0.00854416068 0.4612553 0 0 0.0028000000
## [359,] -0.05274211642 2.0129653 0 0 0.0028000000
## [360,] -0.00377742376 0.6445630 0 0 0.0018666667
## [361,] -0.01158475698 0.6357252 0 0 0.0022666667
## [362,] -0.00008417772 0.4661515 0 0 0.0012000000
## [363,] -0.02499419521 0.9011908 0 0 0.0022666667
## [364,] -0.00443593796 0.3477457 0 0 0.0016000000
## [365,] -0.01049201805 0.5550310 0 0 0.0017333333
## [366,] -0.04586069885 1.4326092 0 0 0.0028000000
## [367,] -0.01109456341 0.5960701 0 0 0.0018666667
## [368,] -0.03920229652 1.3148421 0 0 0.0022666667
## [369,] -0.00129635499 0.5275291 0 0 0.0022666667
## [370,] -0.06081634329 1.8757021 0 0 0.0034666667
## [371,] -0.20830460427 5.8098532 0 0 0.0034666667
## [372,] 0.01657636722 0.9601838 0 0 0.0022666667
## [373,] -0.03899514392 1.7358351 0 0 0.0026666667
## [374,] 0.00365164172 0.6491723 0 0 0.0022666667
## [375,] -0.06984778067 2.0917160 0 0 0.0036000000
## [376,] -0.00452754283 0.7363184 0 0 0.0020000000
## [377,] -0.00537875812 0.3040461 0 0 0.0014666667
## [378,] 0.00344620705 0.2859794 0 0 0.0016000000
## [379,] -0.01373807322 0.5825774 0 0 0.0020000000
## [380,] 0.01395358919 1.1906888 0 0 0.0029333333
## [381,] -0.01840620307 1.0885049 0 0 0.0021333333
## [382,] 0.27126018601 8.1865463 0 0 0.0036000000
## [383,] -0.00070503954 1.2138388 0 0 0.0028000000
## [384,] 0.01373308238 1.0715478 0 0 0.0020000000
## [385,] -0.00267858361 0.7801755 0 0 0.0026666667
## [386,] -0.01676222307 0.6943828 0 0 0.0024000000
## [387,] -0.53665600360 12.2999086 0 0 0.0045333333
## [388,] -0.02291649049 0.9412741 0 0 0.0024000000
## [389,] 0.02127266261 0.9879992 0 0 0.0017333333
## [390,] -0.00346625202 0.6666998 0 0 0.0022666667
## [391,] 0.00132013458 1.6238668 0 0 0.0028000000
## [392,] 0.01196266853 0.8383648 0 0 0.0017333333
## [393,] -0.05426469693 1.3924142 0 0 0.0041333333
## [394,] -0.01865640157 0.5815308 0 0 0.0021333333
## [395,] -0.01084343685 0.6510940 0 0 0.0024000000
## [396,] 0.00323345240 0.5676648 0 0 0.0017333333
## [397,] -0.01327545731 0.7330531 0 0 0.0021333333
## [398,] 0.21992452966 9.1761354 0 0 0.0033333333
## [399,] 0.07591143163 5.3299688 0 0 0.0028000000
## [400,] -0.01186001476 0.9673010 0 0 0.0020000000
## [401,] -0.00361229405 0.6028003 0 0 0.0024000000
## [402,] -0.03123105574 1.0227672 0 0 0.0022666667
## [403,] -0.20809723888 6.9799452 0 0 0.0037333333
## [404,] -0.03460902861 1.3838552 0 0 0.0024000000
## [405,] 0.02214186052 1.1211032 0 0 0.0028000000
## [406,] -0.00603835458 0.5421071 0 0 0.0028000000
## [407,] -0.00324938157 0.3575859 0 0 0.0016000000
## [408,] -0.00626684148 0.8303860 0 0 0.0029333333
## [409,] -0.06835406643 1.5391718 0 0 0.0038666667
## [410,] -0.00690732290 0.3817615 0 0 0.0016000000
## [411,] -0.01980612285 0.9578651 0 0 0.0024000000
## [412,] -0.00673767907 0.5213549 0 0 0.0009333333
## [413,] 0.01489563367 0.8773850 0 0 0.0026666667
## [414,] -0.00180245024 0.4714237 0 0 0.0022666667
## [415,] -0.03341480599 1.0769096 0 0 0.0026666667
## [416,] -0.00326013021 0.6229288 0 0 0.0032000000
## [417,] -0.01017460564 0.9479096 0 0 0.0033333333
## [418,] -0.01942328620 0.7012808 0 0 0.0022666667
## [419,] 0.01150827917 1.5981117 0 0 0.0021333333
## [420,] -0.00872144185 0.7176251 0 0 0.0020000000
## [421,] -0.02406431606 0.9412862 0 0 0.0024000000
## [422,] -0.01574707763 0.9365421 0 0 0.0026666667
## [423,] -0.01555830999 1.7705418 0 0 0.0025333333
## [424,] -0.00005970490 0.6676579 0 0 0.0024000000
## [425,] -0.00230019706 0.4073801 0 0 0.0014666667
## [426,] -0.00147609517 0.2752254 0 0 0.0022666667
## [427,] -0.01950261775 0.8638737 0 0 0.0026666667
## [428,] -0.01425060103 0.6446510 0 0 0.0025333333
## [429,] -0.00705720395 0.5807070 0 0 0.0026666667
## [430,] -0.02623462985 0.8793396 0 0 0.0028000000
## [431,] -0.07339784456 1.9810074 0 0 0.0033333333
## [432,] 0.03737303670 1.3511437 0 0 0.0029333333
## [433,] -0.02254038964 0.9771077 0 0 0.0032000000
## [434,] -0.10306606179 2.8014589 0 0 0.0049333333
## [435,] -0.00839785496 1.3698286 0 0 0.0025333333
## [436,] -0.04891872454 1.3471521 0 0 0.0030666667
## [437,] -0.01298447161 0.8422718 0 0 0.0030666667
## [438,] -0.00197269586 0.3240561 0 0 0.0012000000
## [439,] -0.04287659738 1.5395138 0 0 0.0033333333
## [440,] -0.05220001233 1.6482485 0 0 0.0033333333
## [441,] -0.00850247362 0.5605286 0 0 0.0021333333
## [442,] -0.00638144708 0.5194437 0 0 0.0025333333
## [443,] -0.12408356856 2.7675065 0 0 0.0042666667
## [444,] -0.00905194322 0.7470167 0 0 0.0029333333
## [445,] -0.00812968857 0.7944392 0 0 0.0022666667
## [446,] 0.00001768554 1.0378946 0 0 0.0028000000
## [447,] -0.02218304663 0.9107648 0 0 0.0028000000
## [448,] 0.00092772833 0.6332493 0 0 0.0021333333
## [449,] -0.02568873915 0.9461023 0 0 0.0030666667
## [450,] 0.01392776984 1.4068940 0 0 0.0020000000
## [451,] -0.06793815607 2.2470119 0 0 0.0032000000
## [452,] 0.00436084773 0.7823230 0 0 0.0032000000
## [453,] 0.00366257960 0.4327661 0 0 0.0016000000
## [454,] 0.00268852001 0.4048095 0 0 0.0021333333
## [455,] -0.00480479394 1.3534395 0 0 0.0029333333
## [456,] -0.02585450874 1.0825693 0 0 0.0020000000
## [457,] -0.00902950440 0.4714410 0 0 0.0021333333
## [458,] -0.01651540909 0.6028711 0 0 0.0020000000
## [459,] -0.00387793363 0.3980457 0 0 0.0021333333
## [460,] -0.00002073313 0.3378363 0 0 0.0021333333
## [461,] 0.00499728073 1.1117109 0 0 0.0032000000
## [462,] 0.00794450185 1.0856762 0 0 0.0028000000
## [463,] -0.02138763373 1.2370726 0 0 0.0032000000
## [464,] 0.01974010483 0.9974529 0 0 0.0026666667
## [465,] 0.00483055702 0.4848990 0 0 0.0021333333
## [466,] -0.00646321734 0.5308886 0 0 0.0021333333
## [467,] -0.10459797360 4.1466932 0 0 0.0033333333
## [468,] 0.01493883877 1.0473745 0 0 0.0021333333
## [469,] 0.02249751131 0.9571220 0 0 0.0025333333
## [470,] -0.00837668430 0.4257489 0 0 0.0017333333
## [471,] -0.03773395461 1.0908663 0 0 0.0033333333
## [472,] 0.01837954520 0.8211637 0 0 0.0022666667
## [473,] -0.04562452911 1.0408843 0 0 0.0030666667
## [474,] -0.01371617775 0.8715611 0 0 0.0022666667
## [475,] -0.06018091995 1.5702828 0 0 0.00373333334.2.3 Predict fitted values for each individual
pred.npb <- predict(fit.npb)
fittedvals <- pred.npb$fitted.vals4.2.4 Plot predicted outcomes against “measured” outcomes
plot(fittedvals, Y)
abline(a = 0, b = 1, col = "red")